Figures
NTJ
5/14/2019
Pp.Calc = function(SCA.Sub,Rep) { # Calculates segment-specific Pp values
# Argument = a matrix with rows of cells and columns (PB2...NS,var.PB2,var.PB1,WT,var)
# Need 3 things:
# Number var/Helper+ = A
# Number WT/Test+ = B
# Segment-specific segment counts = Ci
# MOI = B/A
# WT segments can only be propagated in var+ cells, so we consider them the denominator.
# Pp[i] = ln(1 - C[i]/A) / ln(1 - B/A)
require(Matrix)
require(dplyr)
PB2.Index = which(colnames(SCA.Sub)=="PB2")
# Extract Delay = 0 (data where Pan/99-WT and Pan/99-var viruses were added simultaneously) and analyze
Delay0 = SCA.Sub %>% filter(Delay == 0)
A = sum(Delay0[,"var"]) # A = WT+ Cells
B = sum(Delay0[,"WT"]) # B = var+ Cells
MOI = B / A # Here, "MOI"" refers to the fraction of WT+ cells
C = colSums(Delay0[,PB2.Index:(PB2.Index + 7)]) # Counts of each genome segment from WT virus
PP.Sub=log(1 - C / A) / log(1 - MOI)
PP.Sub["Avg"] = prod(PP.Sub) ^ (1 / 8) # Use geometric mean, since product is important
PP.Sub["WT.Segments"] = sum(C)
PP.Sub["P3.NP"] = sum(Delay0[,"P3.NP"]) / B # Fraction of cells that contain WT PB2, PB1, PA, and NP segments
PP.Sub["Segments.Per.Cell"] = sum(C) / (MOI * A) # Average Number of Segments / Cell
PP.Sub["Delay"] = 0 # No delay between adding WT virus and var virus
PP.Sub["Initial.MOI"] = MOI # In the event of a delay between WT-var, the fraction WT+ in the simultaneous co-inoculation
PP.Sub["Time.MOI"] = B / A # In the event of a delay between WT-var, the fraction WT+ at the delay
PP.Sub["Expt"] = SCA.Sub[1,"Expt"] # Experimental Replicate (Name)
PP.Sub["Rep"] = Rep # Experimental Replicate (Number)
PP.Sub
}
Pp.Calc.2 = function(Test) { # Calculates segment-specific Pp values
# Argument = a matrix with rows of cells and columns (PB2...NS,var.PB2,var.PB1,WT,var)
# Need 3 things:
# Number var/Helper+ = A
# Number WT/Test+ = B
# Segment-specific segment counts = Ci
# MOI = B/A
# WT segments can only be propagated in var+ cells, so we consider them the denominator.
# Pp[i] = ln(1 - C[i]/A) / ln(1 - B/A)
require(Matrix)
PB2.Index = which(colnames(Test)=="PB2")
A=sum(Test[,"var"])
B=sum(Test[,"WT"])
C=colSums(Test[,PB2.Index:(PB2.Index + 7)])
PP.Test=log(1-C/A) / log(1-B/A)
PP.Test["Avg"]=prod(PP.Test)^(1/8) # Use geometric mean, since product is important
PP.Test["MOI"]=B/A
PP.Test
}
Virion.By.Genotype.Prob=function(SCA.Sub,Virion.Target) {
vt=Virion.Target # Function is designed to calculate p(1 virion | Genotype, but could be used for an arbitrary number of virions)
PP=Pp.Calc.2(SCA.Sub) # Calculate Pp values
MOI=PP["MOI"]
PP=PP[1:8] # Ignore the "Avg" and "MOI" parts
V.Num=20 #
Cond.Probs=matrix(nrow=256,
ncol=V.Num,
data=0)
Gens=lapply(X=1:255,FUN=IntToGenotype)
Gens=matrix(unlist(Gens,use.names=FALSE),ncol=8,byrow=TRUE)
SF.Probs=matrix(nrow=2,
ncol=8,
data=0)
for (v in 1:V.Num) {
# pgeom (prob = Pp,q = v - 1) = Cumulative Prob of Success after v trials
SF.Probs[2,] = pgeom(prob = PP,q = v-1) #Success Probabilities for each segment (i.e., if a cell is infected with v-1 virions, how likely is it that at least one copy of a given segment is delivered?)
SF.Probs[1,] = 1 - SF.Probs[2,] #Failue Probabilities for each segment (i.e., if a cell is infected with v-1 virions, how likely is it that each virion fails to deliver the segment in question?)
#Each gene constellation or "genotype" represents a distinct combination of WT segments being present and absent. For 8 segments, 2^8 = 256 different genotypes.
Gen.Probs=rep(0,256)
for (gen in 0:255) {
GT = IntToGenotype(gen) # For example, 1 = 00000001, 42 = 00101010, 255 = 11111111
call = cbind(GT+1,1:8)
Gen.Probs[gen + 1] = prod(SF.Probs[call])
}
# Cond.Probs = p(Genotype | # Virions)
Cond.Probs[,v]=Gen.Probs
}
Cond.Probs=cbind(c(1,rep(0,255)),Cond.Probs)
# Cond.Probs = p(Genotype | Virion #) p(B|A)
# Virion.Probs = p(Virion #) p(A)
# Genotype.Probs = p(Genotype) p(B)
# Bayes Theorem: p(A|B) = p(B|A) * p(A) / p(B)
# p(1 virion | Genotype) = p(Genotype | 1 virion) * p(1 virion) / p(Genotype)
# "MOI" calculated in Pp Analysis is actually p (X >= 1)
# Need to calculate lambda to parameterize Poisson distribution
lambda=-log(1 - MOI)
# Need the probability of receiving 0, 1, 2, 3, 4... virions
Virion.Probs=matrix(data=rep(dpois(lambda=lambda,x=0:V.Num),each=256),
ncol=V.Num+1,
nrow=256)
# Genotype.Probs = Sum(i = 0 - 20) (p(Genotype | i Virions) * p(i Virions))
Genotype.Probs = rowSums(Cond.Probs * Virion.Probs)
# p(1 virion | Genotype) = p(Genotype | 1 virion) * p(1 virion) / p(Genotype)
# vt + 1 is an offset, because column 1 contains p(0 virions), column 2 contains p(1 virion), etc.
V.Probs=Cond.Probs[,vt+1] * Virion.Probs[,vt+1] / Genotype.Probs
V.Probs
}
GenotypeToInt = function(Genotype) {
sum(Genotype*2^(0:7))
# For example, 1 = 00000001 (Only NS present)
# 42 = 00101010 (PA, NA, M segments present)
# 255 = 11111111 (all segments present)
}
IntToGenotype = function(Int) {
as.numeric(intToBits(Int)[1:8])
# For example, 1 = 00000001 (Only NS present)
# 42 = 00101010 (PA, NA, M segments present)
# 255 = 11111111 (all segments present)
}
SCA.Cov = function(Test.SCA) { #
PB2.Index = which(colnames(Test.SCA)=="PB2")
Flu.Genes = colnames(Test.SCA)[PB2.Index:(PB2.Index+7)]
Expts = unique(Test.SCA$Expt)
Test.SCA$Weight=0 # To be calculated later, rescaled p1.Virion (sums to 1 for each experiment)
Test.SCA$p1.Virion=0 #To be calculated later
# Calculate Genotypes of each cell. 0 = NO WT segments, 255 = All 8 WT segments
Test.SCA$Genotype=apply(X=Test.SCA[,PB2.Index:(PB2.Index+7)],MARGIN=1,FUN=GenotypeToInt)
Pp.Mat=matrix(data = 0,
nrow = length(Expts),
ncol = 10)
colnames(Pp.Mat) = c(Flu.Genes,"Avg","MOI")
rownames(Pp.Mat) = str_c("Exp_",1:length(Expts))
for (i in 1:length(Expts)) {
# Calculate Pp
Pp.Mat[i,]=Pp.Calc.2(Test = subset(Test.SCA,Expt==Expts[i]))
# Calculate p(1 Virion | Genotype), and Weight (relative p(1|Genotype), sums to 1 in each Expt)
Single.Hit <- Virion.By.Genotype.Prob(subset(Test.SCA,Expt==Expts[i]),
Virion.Target = 1)
Test.SCA[Test.SCA[,"Expt"]==Expts[i],"p1.Virion"] <- Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1]
Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Weight"] <- Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1] / sum(Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1])
}
# After calculating p(1 virion) for each cell, use the cov.wt function to calculate the pairwise correlation between each pair of segments.
SCA.Ass = cov.wt(x = Test.SCA[,(PB2.Index:(PB2.Index + 7))], wt = Test.SCA[,"p1.Virion"],
cor = TRUE, center = TRUE)
return(SCA.Ass)
}
Weight.Calc = function(Test.SCA) {
# Calculate p(1 virion) for each cell, to be used as a weighting factor in the pairwise association analysis
PB2.Index = which(colnames(Test.SCA)=="PB2")
Flu.Genes = colnames(Test.SCA)[PB2.Index:(PB2.Index+7)]
Expts = unique(Test.SCA$Expt)
Test.SCA$Weight=0 # To be calculated later, rescaled p1.Virion (sums to 1 for each experiment)
Test.SCA$p1.Virion=0 #To be calculated later
# Calculate Genotypes of each cell. 0 = NO WT segments, 255 = All 8 WT segments
Test.SCA$Genotype=apply(X=Test.SCA[,PB2.Index:(PB2.Index+7)],MARGIN=1,FUN=GenotypeToInt)
Pp.Mat=matrix(data = 0,
nrow = length(Expts),
ncol = 10)
colnames(Pp.Mat) = c(Flu.Genes,"Avg","MOI")
rownames(Pp.Mat) = str_c("Exp_",1:length(Expts))
for (i in 1:length(Expts)) {
# Calculate Pp
Pp.Mat[i,]=Pp.Calc.2(Test = subset(Test.SCA,Expt==Expts[i]))
# Calculate p(1 Virion | Genotype), and Weight (relative p(1|Genotype), sums to 1 in each Expt)
Single.Hit <- Virion.By.Genotype.Prob(subset(Test.SCA,Expt==Expts[i]),
Virion.Target = 1)
Test.SCA[Test.SCA[,"Expt"]==Expts[i],"p1.Virion"] <- Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1]
Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Weight"] <- Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1] / sum(Single.Hit[Test.SCA[Test.SCA[,"Expt"]==Expts[i],"Genotype"]+1])
}
return(Test.SCA)
}
SCA.Data = read.csv(file.path(Data.Path,"SCA_Full_Data.csv")) %>%
filter(Cell.Type =="Test") %>%
mutate(WT.Segments = PB2 + PB1 + PA + HA + NP + NA. + M + NS,
WT = sign(WT.Segments),
var = sign(Help.PB2 * Help.PB1),
P3.NP = sign(PB2 * PB1 * PA * NP))
## Calculate Pp values algorithmically ----
SCA.Expts = unique(SCA.Data$Expt)
Pp.Summary = SCA.Data %>%
filter(Expt == SCA.Expts[1]) %>%
Pp.Calc(Rep = 1)
for (i in 2:length(SCA.Expts)) {
Pp.Summary = rbind(Pp.Summary,
SCA.Data %>%
filter(Expt == SCA.Expts[i]) %>%
Pp.Calc(Rep = i))
}
rownames(Pp.Summary) = NULL
Pp.Summary = as.data.frame(Pp.Summary)
write.csv(Pp.Summary,file = file.path(Data.Path,"Pp_Summary.csv"))
# Pp vs. Length (0 hrs) ----
Rep.Num = length(SCA.Expts)
Exp.Pp = Pp.Summary %>% filter(Delay == 0) %>% select(PB2:NS,Rep)
colnames(Exp.Pp)[which(colnames(Exp.Pp)=="NA.")]="NA"
Exp.Pp=melt(Exp.Pp,id=c("Rep"))
Exp.Pp=Exp.Pp[1:(Rep.Num * 8),]
Segment.Names=c("PB2","PB1","PA","HA","NP","NA","M","NS")
Segment.Lengths=c(2341,2341,2233,1778,1565,1413,1027,890)
Lengths=rep(c(2341,2341,2233,1778,1565,1413,1027,890),each=length(unique(Exp.Pp$Rep)))
Exp.Pp=cbind(Exp.Pp,Lengths)
colnames(Exp.Pp)[which(colnames(Exp.Pp)=="value")]="Pp"
colnames(Exp.Pp)[which(colnames(Exp.Pp)=="variable")]="Segment"
Exp.Pp$Pp=as.numeric(Exp.Pp$Pp)
Exp.Pp$Lengths=as.numeric(Exp.Pp$Lengths)
Exp.Pp$Rep = factor(Exp.Pp$Rep,levels = 1:Rep.Num)
# Segment Distribution ----
# Plot histogram of segments/cell (only WT+ cells) and ECDF
x = SCA.Data %>%
mutate(WT.Segments = as.numeric(WT.Segments)) %>%
filter(WT == 1,
Delay == 0) %>%
dplyr::select(WT.Segments)
WT_at_least = rep(0,8)
WT_cdf = rep(0,8)
WT_hist = table(x) / nrow(x)
Max_Freq = max(WT_hist)
for (i in 1:8) {
WT_at_least[i] = sum(x >= i) / nrow(x) * Max_Freq
WT_cdf[i] = sum(x <= i) / nrow(x) * Max_Freq
}
Segment_Dist = cbind(1:8,WT_at_least,WT_cdf,WT_hist) %>% data.frame
colnames(Segment_Dist) = c("Segments","At_Least","ECDF","Freq")
ggplot() +
geom_bar(data = Segment_Dist,
aes(x = Segments,
y = Freq),
stat = "identity",
fill = "lightgray",
color = "black") +
geom_line(data = Segment_Dist,
aes(x = Segments,
y = ECDF),
color = "red",
lwd = 1.2) +
labs(x = "Segments / Cell",
y = "Individual Frequency") +
scale_x_continuous(breaks = 1:8) +
scale_y_continuous(sec.axis = sec_axis(~. / Max_Freq, name = "Cumulative Frequency")) +
theme(text=element_text(size=12,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.title.y.right = element_text(color = "red"),
axis.text.y.right = element_text(color = "red"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave("Misc_Plots/Segment_Distribution.pdf",
dpi = 300,
width = 5,
height = 4,
unit = "in")
# Stock Plot ----
Fig.1 = ggplot() +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = c(2,2))
# Figure 1A: Pp by Segment ----
Pp_Error = Exp.Pp %>%
group_by(Segment) %>%
dplyr::summarise(Pp_Mean = mean(Pp),
Pp_SE = sd(Pp) / sqrt(length(Pp)))
Fig.1 +
geom_jitter(data = Exp.Pp,
aes(x = Segment,
y = Pp,
color = Rep),
size = 2) +
geom_crossbar(data = Pp_Error,
aes(x = Segment,
y = Pp_Mean,
ymin = Pp_Mean - 1.96 * Pp_SE,
ymax = Pp_Mean + 1.96 * Pp_SE),
fill = "darkgray",
alpha = 0.5) +
scale_y_continuous(limits=c(0,1.05),breaks=0:5/5) +
annotate("text",x = 3, y = 0.99,size = 5,label = TeX("\\textbf{$P_P$ = 0.58 (0.57 - 0.59)}")) +
labs(y=expression(bold(paste("Probability of Detection (P"[bold("P")],")"))),
x="Segment",
color="Replicate") +
scale_color_tableau(palette = "Tableau 20")## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
ggsave("Figures/1A_Pp_By_Segment.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'write.csv(Exp.Pp %>% dplyr::select(Rep,Segment,Pp), file = file.path(Data.Path,"Fig1A_Data.csv"), row.names = FALSE)
# Legend
Legend.1AB = Fig.1 +
geom_jitter(data = Exp.Pp,
aes(x = Segment,
y = Pp,
color = Rep),
size = 2) +
geom_crossbar(data = Pp_Error,
aes(x = Segment,
y = Pp_Mean,
ymin = Pp_Mean - 1.96 * Pp_SE,
ymax = Pp_Mean + 1.96 * Pp_SE),
fill = "darkgray",
alpha = 0.5) +
scale_color_tableau(palette = "Tableau 20") +
scale_y_continuous(limits=c(0,1.05),breaks=0:5/5) +
annotate("text",x = 4, y = 0.99,size = 5,label = TeX("\\textbf{$P_P$ = 0.58 (0.57 - 0.59)}")) +
labs(y=expression(bold(paste("Probability of Detection (P"[bold("P")],")"))),
x="Segment",
color="Replicate") +
theme(legend.position = "right")
ggsave(g_legend(Legend.1AB),
file = "Figures/1AB_Legend.pdf",
dpi = 300,
width = 2,
height = 5,
unit = "in")## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'# Predicted Reassortment based on Pp values----
FM.Data = read.table(file = file.path(Data.Path,"Fonville_Marshall_Reassortment_Data.csv"),header=TRUE,sep=",")
# Predictions generated by experimentally measured Pp values
Pp_Reassortment_Predict = read.csv(file = file.path(Data.Path,"Exp_Pp_Reassortment.csv"),header = TRUE, sep = ",")
# Expected reassortment frequencies if genomes were always complete (Pp = 1)
Pp1_Reassortment_Predict = read.csv(file = file.path(Data.Path,"Pp1_Reassortment.csv"),header = TRUE, sep = ",")
X="Expressing.HA"
Y="Reassortant.Percent"
Fig.1 +
geom_line(data=Pp_Reassortment_Predict,aes(y=Reassortant.Percent,x=Expressing.HA,group=as.factor(Rep),color=as.factor(Rep)),lwd=1.2) +
geom_line(data = Pp1_Reassortment_Predict, aes(x = Expressing.HA, y = Reassortant.Percent), lty = 2, lwd = 1.2) +
geom_point(data=FM.Data,aes_string(x=X,y=Y),size=4,col="black") +
annotate("text",x = 75, y = 15,size = 6,label = TeX("\\textbf{$\\P_P$ = 1}")) +
labs(x=expression(bold("% Cells HA"^"+")),
y="% Reassortment",
color="Rep") +
scale_color_tableau(palette = "Tableau 20") +
guides(color = FALSE) +
scale_x_continuous(limits=c(0,100)) +
scale_y_continuous(limits=c(0,100))## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
ggsave("Figures/1B_Predicted_Reassortment.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'write.csv(rbind(Pp_Reassortment_Predict %>% dplyr::select(Rep,Expressing.HA,Reassortant.Percent),
Pp1_Reassortment_Predict %>% mutate(Rep = "Pp_1") %>% dplyr::select(Rep,Expressing.HA,Reassortant.Percent),
FM.Data %>% mutate(Rep = "Fonville_Marshall_Experiments")), file = file.path(Data.Path,"Fig1B_Data.csv"), row.names = FALSE)
# Figure 1C: Pairwise Segment Correlations ----
# Calculate Genotype #s and p(1 virion | Genotype)
SCA.Weights = Weight.Calc(Test.SCA = SCA.Data %>% filter(Delay == 0,
Cell.Type =="Test",
var == 1))
SCA.Ass = SCA.Cov(Test.SCA = SCA.Data %>% filter(Delay == 0,
Cell.Type =="Test",
var == 1))
# Plotting Weights ----
ggplot() +
geom_jitter(data = SCA.Weights,
aes(x = WT.Segments,
y = p1.Virion))
ggplot() +
geom_density_ridges(data = SCA.Weights %>% filter(WT.Segments >= 1),
aes(x = p1.Virion,
y = factor(WT.Segments,levels = rev(1:8))),
fill = "navy",
alpha = 0.8,
scale = 2) +
scale_x_continuous(limits = c(-0.05,1.08),
breaks = c(0,0.25,0.5,0.75,1.0)) +
labs(x = "p (1 Virion)",
y = "# WT Segments Detected") +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())## Picking joint bandwidth of 0.0273
ggsave("Figures/Supp1_Weight_by_Segment.pdf",
dpi = 300,
width = 5,
height = 4,
unit = "in")## Picking joint bandwidth of 0.0273write.csv(SCA.Weights %>% filter(WT.Segments >= 1) %>% select(Expt,Cell.ID:Help.PB1,p1.Virion,Genotype), file = file.path(Data.Path,"Supp1_Data.csv"), row.names = FALSE)
# Plotting Correlation ----
Flu.Genes = c("PB2","PB1","PA","HA","NP","NA","M","NS")
#Flu.Genes[6]="NA"
Sum.Names=c("Index","Pair","X.Num","Y.Num","X","Y","r","r2","p")
Sum.df=matrix(data=0,nrow=64,ncol=length(Sum.Names))
colnames(Sum.df)=Sum.Names
Sum.df[,"Index"]=1:64
Sum.df[,"X"]=rep(Flu.Genes,each=8)
Sum.df[,"Y"]=rep((Flu.Genes),times=8)
Sum.df[,"r"]=as.numeric(SCA.Ass$cor)
Sum.df=as.data.frame(Sum.df)
Sum.df$r = as.numeric(SCA.Ass$cor)
Sum.df$r2 = as.numeric(Sum.df$r * Sum.df$r)
Sum.df$X=factor(Sum.df$X,levels=Flu.Genes)
Sum.df$Y=factor(Sum.df$Y,levels=rev(Flu.Genes))
Sum.df$p = pmin((1 - pchisq(Sum.df$r2 * sum(SCA.Weights$p1.Virion), df = 3)) * 28, 1)
for (row in 1:8) {
for (col in 1:8) {
Pair.Index=(row-1)*8+col
if (col<=row) {
Sum.df[Pair.Index,"r2"]=-1
}
}
}
Sum.df=Sum.df[Sum.df[,"X"]!="NS",]
Sum.df=Sum.df[Sum.df[,"Y"]!="PB2",]
Sum.df=Sum.df[Sum.df[,"r2"]!=-1,]
ggplot() +
coord_equal() +
geom_tile(data=Sum.df,aes(x=X,
y=(Y),
fill=r)) +
geom_text(data = Sum.df %>% filter(p < 0.05),
aes(x = X,
y = Y,
label = round(r,2)),
color = "yellow",
size = 5) +
scale_fill_viridis(end = max(Sum.df$r), option = "inferno") +
labs(x=NULL,
y=NULL,
fill=expression(bold("r")^bold(""))) +
theme(legend.position="NA",
text=element_text(size=14,face="bold"),
plot.title=element_text(size=rel(2),face="bold"),
legend.title.align=3,
legend.title=element_text(size=rel(1.5),face="bold",vjust=0),
legend.text=element_text(size=rel(1.2),face="bold",vjust=0),
legend.key.width=unit(2,"cm"),
axis.text.x = element_text(size=rel(1.5),vjust = .5,hjust=.5,face="bold",color = "black", margin=margin(5,0,0,0), angle=0),
axis.text.y = element_text(size=rel(1.5),face="bold",color = "black"),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line=element_blank(),
strip.text.x = element_text(size=rel(1),angle=0,face="bold"),
strip.text.y = element_text(size=rel(1),angle=0,face="bold"),
strip.background = element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background = element_blank())
ggplot() +
coord_equal() +
geom_tile(data=Sum.df,aes(x=X,
y=(Y),
fill=r2)) +
geom_text(data = Sum.df %>% filter(p < 0.05),
aes(x = X,
y = Y,
label = round(r2,2)),
color = "yellow",
size = 5) +
scale_fill_viridis(end = max(Sum.df$r2), option = "inferno") +
labs(x=NULL,
y=NULL,
fill=expression(bold("r")^bold(""))) +
theme(legend.position="NA",
text=element_text(size=14,face="bold"),
plot.title=element_text(size=rel(2),face="bold"),
legend.title.align=3,
legend.title=element_text(size=rel(1.5),face="bold",vjust=0),
legend.text=element_text(size=rel(1.2),face="bold",vjust=0),
legend.key.width=unit(2,"cm"),
axis.text.x = element_text(size=rel(1.5),vjust = .5,hjust=.5,face="bold",color = "black", margin=margin(5,0,0,0), angle=0),
axis.text.y = element_text(size=rel(1.5),face="bold",color = "black"),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line=element_blank(),
strip.text.x = element_text(size=rel(1),angle=0,face="bold"),
strip.text.y = element_text(size=rel(1),angle=0,face="bold"),
strip.background = element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background = element_blank())
ggsave("Figures/1C_Pairwise_Correlation_r2.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")
write.csv(Sum.df %>% dplyr::select(X,Y,r,r2,p), file = file.path(Data.Path,"Fig1C_Data.csv"), row.names = FALSE)
Legend.1C = ggplot() +
coord_equal() +
geom_tile(data=Sum.df %>% mutate(r2 = r2 / max(r2)),aes(x=X,
y=(Y),
fill=r2)) +
scale_fill_viridis(begin = 0, end = 1, option = "inferno") +
labs(x=NULL,
y=NULL,
fill=expression(bold("r")^bold("2"))) +
theme(text=element_text(size=14,face="bold"),
plot.title=element_text(size=rel(2),face="bold"),
legend.title.align=3,
legend.title=element_text(size=rel(1.5),face="bold",vjust=0),
legend.text=element_text(size=rel(1.2),face="bold",vjust=0),
legend.key.width=unit(2,"cm"),
legend.position = "bottom",
axis.text.x = element_text(size=rel(1.5),vjust = .5,hjust=.5,face="bold",color = "black", margin=margin(5,0,0,0), angle=0),
axis.text.y = element_text(size=rel(1.5),face="bold",color = "black"),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line=element_blank(),
strip.text.x = element_text(size=rel(1),angle=0,face="bold"),
strip.text.y = element_text(size=rel(1),angle=0,face="bold"),
strip.background = element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background = element_blank())
ggsave(g_legend(Legend.1C),
file = "Figures/1C_Legend_r2.pdf",
dpi = 300,
width = 5,
height = 2,
unit = "in")
# Pp CI Calculation ----
Pp.Summary %>% select(PB2:NS) %>%
melt %>%
rename(Pp = value,
Segment = variable) %>%
summarise(Pp.Mean = mean(Pp),
Pp.se = sd(Pp)/length(Pp)) %>%
mutate(Pp.Min = Pp.Mean - 1.96 * Pp.se,
Pp.Max = Pp.Mean + 1.96 * Pp.se)## No id variables; using all as measure variables## Pp.Mean Pp.se Pp.Min Pp.Max
## 1 0.5858705 0.001215225 0.5834886 0.5882523Pp.Summary %>% select(Avg) %>%
summarise(Pp.Mean = mean(Avg),
Pp.se = sd(Avg)/length(Avg)) %>%
mutate(Pp.Min = Pp.Mean - 1.96 * Pp.se,
Pp.Max = Pp.Mean + 1.96 * Pp.se)## Pp.Mean Pp.se Pp.Min Pp.Max
## 1 0.5758761 0.005454031 0.5651862 0.586566# 0.58 (0.57 - 0.59) 2 digitsCellInfection=function(Cell.Num,MOI,Cell.Type,Pp) {
set.seed(627) #
Cell.Names = c(
"MOI","Cell.Type",
"Cell","Infected","Co.Infected","Productive","Productive.Co.Infected","Poly","Surface.HA","Inf.A","Inf.B",
"pParent.A","pParent.B","pReassort","pHN.Reassort",
"PB2","PB1","PA","HA","NP","NA","M","NS",
"A.PB2","A.PB1","A.PA","A.HA","A.NP","A.NA","A.M","A.NS",
"B.PB2","B.PB1","B.PA","B.HA","B.NP","B.NA","B.M","B.NS",
"Copy.A.PB2","Copy.A.PB1","Copy.A.PA","Copy.A.HA","Copy.A.NP","Copy.A.NA","Copy.A.M","Copy.A.NS",
"Copy.B.PB2","Copy.B.PB1","Copy.B.PA","Copy.B.HA","Copy.B.NP","Copy.B.NA","Copy.B.M","Copy.B.NS",
"Copy.PB2","Copy.PB1","Copy.PA","Copy.HA","Copy.NP","Copy.NA","Copy.M","Copy.NS",
"Six","HA.Neg","NS.Neg"
)
Cells = matrix(nrow = Cell.Num,
ncol = length(Cell.Names),
data = 0)
colnames(Cells) = Cell.Names
Cells[, "Cell"] = seq(1, Cell.Num)
Cells[, "MOI"] = MOI
withA = sort(ceiling(runif(n = MOI * Cell.Num) * Cell.Num))
withB = sort(ceiling(runif(n = MOI * Cell.Num) * Cell.Num))
x = table(withA)
x_index = as.numeric(names(x))
x_values = as.numeric(x)
Cells[x_index, "Inf.A"] = x_values
x = table(withB)
x_index = as.numeric(names(x))
x_values = as.numeric(x)
Cells[x_index, "Inf.B"] = x_values
Cells[, "Co.Infected"] = (Cells[, "Inf.A"] * Cells[, "Inf.B"] > 0)
Cells[, "Infected"] = (Cells[, "Inf.A"] + Cells[, "Inf.B"] > 0)
Cells[, "pParent.A"] = 1
Cells[, "pParent.B"] = 1
Cells[, "pHN.Reassort"] = 1
A.PB2.Index = which((colnames(Cells)) == "A.PB2")
PB2.Index = which((colnames(Cells)) == "PB2")
for (segment in 1:8) {
x = table(sort(sample(
withA,
size = length(withA) * Pp[segment],
replace = FALSE
)))
x_index = as.numeric(names(x))
x_values = as.numeric(x)
Cells[x_index, A.PB2.Index + segment - 1] = x_values > 0 #Boolean Presence/Absence of Segment, A strain
Cells[x_index, A.PB2.Index + segment - 1 + 16] = x_values #Copy Number of each segment, A strain
x = table(sort(sample(
withB,
size = length(withB) * Pp[segment],
replace = FALSE
)))
x_index = as.numeric(names(x))
x_values = as.numeric(x)
Cells[x_index, A.PB2.Index + segment - 1 + 8] = x_values > 0 #Boolean Presence/Absence of Segment, B strain
Cells[x_index, A.PB2.Index + segment - 1 + 16 + 8] = x_values #Copy Number of each segment, B strain
Cells[, A.PB2.Index + segment - 1 + 16 + 8 + 8] = Cells[, A.PB2.Index + segment - 1 + 16] + Cells[, A.PB2.Index + segment - 1 + 16 + 8] #Copy number of each segment, irrespective of source
Cells[, PB2.Index + segment - 1] = (Cells[, A.PB2.Index + segment - 1] +
Cells[, A.PB2.Index + segment - 1 + 8]) > 0 #Presence/Absence of Gene, irrespective of source
Cells[, "pParent.A"] = Cells[, "pParent.A"] * (Cells[, A.PB2.Index + segment -
1 + 16] / Cells[, A.PB2.Index + segment - 1 + 16 + 8 + 8])
Cells[, "pParent.B"] = Cells[, "pParent.B"] * (Cells[, A.PB2.Index + segment -
1 + 16 + 8] / Cells[, A.PB2.Index + segment - 1 + 16 + 8 + 8])
}
Cells[, "Poly"] = Cells[, "PB2"] * Cells[, "PB1"] * Cells[, "PA"] * Cells[, "NP"]
Cells[, "Surface.HA"] = Cells[, "Poly"] * Cells[, "HA"]
Cells[, "Productive"] = Cells[, "Surface.HA"] * Cells[, "M"] * Cells[, "NS"] *
Cells[, "NA"]
Cells[, "Productive.Co.Infected"] = Cells[, "Productive"] * Cells[, "Co.Infected"]
Cells[, "pReassort"] = 1 -
Cells[, "pParent.A"] -
Cells[, "pParent.B"]
Cells[, "pHN.Reassort"] = 1 -
(Cells[, "Copy.A.HA"] / Cells[, "Copy.HA"]) * (Cells[, "Copy.A.NA"] / Cells[, "Copy.NA"]) -
(Cells[, "Copy.B.HA"] / Cells[, "Copy.HA"]) * (Cells[, "Copy.B.NA"] / Cells[, "Copy.NA"])
Cells
}
### Parameters ----
Exp.Pp.2 = read.csv(file.path(Data.Path,"Pp_Summary.csv"),na.strings=c("","-")) %>% filter(Delay == 0)
PB2.Index = which(colnames(Exp.Pp.2) == "PB2")
Pp = Exp.Pp.2[,PB2.Index:(PB2.Index + 7)]
Rep.Num=nrow(Pp)
rownames(Pp)=as.character(1:Rep.Num)
MOI=matrix(data=rep(c(.001,.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1,1.5,2,5,10),Rep.Num), #Range of MOIs
byrow=TRUE,nrow=Rep.Num)
rownames(MOI)=rownames(Pp)
Flu.Segments=c("PB2","PB1","PA","HA","NP","NA","M","NS")
MOI.Num=ncol(MOI)
Cell.Num = 100000 # Usually inoculate 1,000,000 cells, but similar results obtained by simulating 100,000 cells
Burst.Size = 986
Rep.Num=nrow(Pp)
### Initialize Summary Data Frame ----
Cell.Summary.Names=c(
# Model Parameters
"Rep","MOI",
# Outcomes
"Productive", "Productive.Co.Infected","Expressing.HA",
"Virion.Number","Parent.A.Percent","Parent.B.Percent","Reassortant.Percent","Reassortant.HN.Percent",
"Parent.A.Number","Parent.B.Number","Reassortant.Number","Reassortant.HN.Number"
)
Cell.Summary=matrix(nrow=(MOI.Num*Rep.Num),
ncol=length(Cell.Summary.Names),
data=0)
colnames(Cell.Summary)=Cell.Summary.Names
Cell.Summary=as.data.frame(Cell.Summary)
Cell.Summary$Rep=rownames(Pp)
Complete.Summary = Cell.Summary[1:MOI.Num,]
Summary.Index = 1
Pp = as.matrix(Pp)
# Simulation ----
Sim.Start = Sys.time()
for (moi in 1:MOI.Num) {
for (cell in 1:Rep.Num) {
Infected.Cells = CellInfection(MOI = MOI[cell,moi],
Cell.Num = Cell.Num,
Pp = Pp[cell,])
Cell.Summary[Summary.Index,"Productive"] = mean(Infected.Cells[,"Productive"]) * 100
Cell.Summary[Summary.Index,"Productive.Co.Infected"] = mean(Infected.Cells[,"Productive.Co.Infected"]) * 100
Cell.Summary[Summary.Index,"Expressing.HA"] = mean(Infected.Cells[,"Surface.HA"]) * 100
Cell.Summary[Summary.Index,"Rep"] = rownames(Pp)[cell]
Cell.Summary[Summary.Index,"MOI"] = MOI[cell,moi]
Cell.Summary[Summary.Index,"Virion.Number"] = Burst.Size * Cell.Summary[Summary.Index,"Productive"]
Infected.Cells=Infected.Cells[Infected.Cells[,"Productive"]==1,,drop=FALSE]
Cell.Summary[Summary.Index,"Parent.A.Percent"] = mean(Infected.Cells[,"pParent.A"]) * 100
Cell.Summary[Summary.Index,"Parent.B.Percent"] = mean(Infected.Cells[,"pParent.B"]) * 100
Cell.Summary[Summary.Index,"Reassortant.Percent"] = mean(Infected.Cells[,"pReassort"]) * 100
Cell.Summary[Summary.Index,"Reassortant.HN.Percent"] = mean(Infected.Cells[,"pHN.Reassort"]) * 100
Cell.Summary[Summary.Index,"Parent.A.Number"] = Cell.Summary[Summary.Index,"Parent.A.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Cell.Summary[Summary.Index,"Parent.B.Number"] = Cell.Summary[Summary.Index,"Parent.B.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Cell.Summary[Summary.Index,"Reassortant.Number"] = Cell.Summary[Summary.Index,"Reassortant.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Cell.Summary[Summary.Index,"Reassortant.HN.Number"] = Cell.Summary[Summary.Index,"Reassortant.HN.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Summary.Index=Summary.Index+1
}
Infected.Cells = CellInfection(MOI = MOI[1,moi],
Cell.Num = Cell.Num,
Pp = rep(1,8))
Complete.Summary[moi,"Productive"] = mean(Infected.Cells[,"Productive"]) * 100
Complete.Summary[moi,"Productive.Co.Infected"] = mean(Infected.Cells[,"Productive.Co.Infected"]) * 100
Complete.Summary[moi,"Expressing.HA"] = mean(Infected.Cells[,"Surface.HA"]) * 100
Complete.Summary[moi,"Rep"] = 1
Complete.Summary[moi,"MOI"] = MOI[cell,moi]
Complete.Summary[moi,"Virion.Number"] = Burst.Size * Cell.Summary[Summary.Index,"Productive"]
Infected.Cells=Infected.Cells[Infected.Cells[,"Productive"]==1,,drop=FALSE]
Complete.Summary[moi,"Parent.A.Percent"] = mean(Infected.Cells[,"pParent.A"]) * 100
Complete.Summary[moi,"Parent.B.Percent"] = mean(Infected.Cells[,"pParent.B"]) * 100
Complete.Summary[moi,"Reassortant.Percent"] = mean(Infected.Cells[,"pReassort"]) * 100
Complete.Summary[moi,"Reassortant.HN.Percent"] = mean(Infected.Cells[,"pHN.Reassort"]) * 100
Complete.Summary[moi,"Parent.A.Number"] = Cell.Summary[Summary.Index,"Parent.A.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Complete.Summary[moi,"Parent.B.Number"] = Cell.Summary[Summary.Index,"Parent.B.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Complete.Summary[moi,"Reassortant.Number"] = Cell.Summary[Summary.Index,"Reassortant.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
Complete.Summary[moi,"Reassortant.HN.Number"] = Cell.Summary[Summary.Index,"Reassortant.HN.Percent"] * Cell.Summary[Summary.Index,"Virion.Number"]
}
Sim.End = Sys.time()
Sim.Time = Sim.End - Sim.Start
# Analysis ----
X="Expressing.HA"
Y="Reassortant.Percent"
FM.Data=read.csv(file = file.path(Data.Path,"Fonville_Marshall_Reassortment_Data.csv"),header=TRUE,sep=",")
Cell.Summary$Rep = as.numeric(Cell.Summary$Rep)
ggplot() +
geom_line(data=Cell.Summary,aes(y=Reassortant.Percent,x=Expressing.HA,group=as.factor(Rep),color=as.factor(Rep)),lwd=1) +
geom_line(data = Complete.Summary, aes(x = Expressing.HA, y = Reassortant.Percent), lty = 2, lwd = 1) +
geom_point(data=FM.Data,aes_string(x=X,y=Y),size=2,col="black") +
annotate("text",x = 75, y = 15,size = 6,label = TeX("\\textbf{$\\P_P$ = 1}")) +
labs(x=expression(bold("% Cells HA"^"+")),
y="% Reassortment",
color="Rep") +
guides(color = FALSE) +
scale_x_continuous(limits=c(0,100)) +
scale_y_continuous(limits=c(0,100)) +
theme(text=element_text(size=12,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave(file.path(Data.Path,"Figures/1B_Predicted_Reassortment.pdf"),
dpi = 300,
width = 4,
height = 4,
unit = "in")Pp = seq(from = 0.01, to = 1, by = 0.01)
tau = matrix(nrow = 1,
ncol = 9,
data = c(1,0,0,0,0,0,0,0,0))
sub.tau = matrix(nrow = 1,
ncol = 8,
data = c(1,0,0,0,0,0,0,0))
trans = matrix(nrow = 9,
ncol = 9,
data = 0)
ident = diag(8)
one.mat = matrix(nrow = 8,
ncol = 1,
data = 1)
Segment.Probs = matrix(ncol = 5,
nrow = 100 * 100 * 9,
data = 0)
Inf.Unit = rep(0, 100)
Inf.Unit.sd = rep(0,100)
colnames(Segment.Probs) = c("Pp","Virion","Segments","Prob","Expected.Segments")
# Order = Pp -> Virion # -> Segments
Segment.Probs[,"Pp"] = rep(seq(from = 0.01, to = 1, by = 0.01), each = 900)
Segment.Probs[,"Virion"] = rep(1:100, each = 9)
Segment.Probs[,"Segments"] = rep(seq(0,8),100)
Segment.Mat = matrix(nrow = 9,
ncol = 1,
data = 0:8)
index = 1
for (pp in 1:100) {
# Define Transition Matrix
for (i in 0:8) {
for (j in i:8) {
trans[i+1,j+1] = dbinom(x = j - i, size = 8 - i, prob = Pp[pp])
}
trans[9,9] = 1
}
for (v in 1:100) {
Prob.Vector = t(tau %*% (trans %^% v))
Segment.Probs[index:(index + 8),"Prob"] = Prob.Vector
Segment.Probs[index:(index + 8),"Expected.Segments"] = t(Prob.Vector) %*% Segment.Mat
index = index + 9
}
# Calculate Infectious Unit
sub.trans = trans[1:8,1:8]
Q = sub.trans
N = solve(ident - Q)
t = N %*% one.mat
tsq = t^2
Inf.Unit[pp] = sub.tau %*% solve(ident - sub.trans) %*% one.mat
Inf.Unit.sd[pp] = sqrt(((2 * N - ident) %*% t - tsq)[1,1])
}
Segment.Probs = as.data.frame(Segment.Probs)
# Productive Matrix ----
Productive = Segment.Probs %>% filter(Segments == 8)
Productive = Productive[,c(1,2,4,5)]
Productive$Convert = Productive$Prob
index = 1
for (pp in 1:100) {
Productive[(index + 1):(index + 99),"Convert"] =
Productive[(index + 1):(index + 99),"Prob"] - Productive[(index):(index + 98),"Prob"]
index = index + 100
}
Geom = matrix(ncol = 3,
nrow = 100,
data = c(Pp,Inf.Unit,Inf.Unit.sd),
byrow = FALSE) %>%
as.data.frame
colnames(Geom) = c("Pp","Inf.Unit","Inf.Unit.sd")
Exp.Pp = read.csv(file = file.path(Data.Path,"Pp_Summary.csv"),na.strings = c("","-")) %>%
filter(Delay == 0) %>%
dplyr::select(Rep,Avg) %>%
mutate(Pp = as.numeric(round(Avg * 100)))
Exp.Pp$Rep = (Exp.Pp$Rep)
Inf.Unit.Sum = Geom[Exp.Pp[1:13,"Pp"],]
Inf.Unit.Sum$Rep = 1:nrow(Inf.Unit.Sum)
Inf.Dose = mean(Inf.Unit.Sum$Inf.Unit) %>% round(1)
Geom$Particle_PFU = (Geom$Pp ^ 8)
Inf.Unit.Sum$Full = (Inf.Unit.Sum$Pp ^ 8) * 100
Geom$Full.sd = sqrt(Geom$Particle_PFU * (1 - Geom$Particle_PFU))
Inf.Unit.Sum$Full %>% mean## [1] 1.772705Inf.Unit.Sum$Inf.Unit %>% mean## [1] 3.728939Fig.2 = ggplot() +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5, color = "black"),
axis.line.y = element_line(size=0.5, color = "black"),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
# Figure 2A: # Segments from 1 Virion ----
Fig.2 +
geom_line(data = Segment.Probs %>% filter(Virion == 1,
floor(Pp*100) %in% c(10,58,90)),
aes(x = Segments,
y = Prob,
color = Pp,
group = Pp),
lwd = 1) +
labs(x = "# Segments",
y = "Frequency") +
scale_color_viridis(begin = 0.1, end = 0.9)
ggsave("Figures/2A_Segments_from_1_Virion.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")
write.csv(Segment.Probs %>% filter(Virion == 1, floor(Pp*100) %in% c(10,58,90)), file = file.path(Data.Path,"Fig2A_Data.csv"), row.names = FALSE)
# Figure 2B: % Fully Infectious Virions ----
Fig.2 +
geom_line(data = Geom,
aes(x = Pp,
y = Particle_PFU * 100),
lwd = 1) +
geom_ribbon(data = Geom,
aes(x = Pp,
ymin = pmax(0,(Particle_PFU - Full.sd * 1.96) * 100),
ymax = pmin(100,(Particle_PFU + Full.sd * 1.96) * 100)),
alpha = 0.3) +
geom_point(data = Inf.Unit.Sum,
aes(x = Pp,
y = Full,
color = factor(Rep)),
size = 2) +
geom_segment(aes(y = Inf.Unit.Sum$Full %>% mean,
yend = Inf.Unit.Sum$Full %>% mean,
# y = 1.22,
# yend = 1.22,
x = 0,
xend = 0.58),
linetype = 2) +
xlab(bquote(bold(P[P]))) +
scale_color_tableau(palette = "Tableau 20") +
scale_y_continuous(breaks = c(1.2,25,50,75,100)) +
# annotate("text",x = 0.3, y = 80,size = 5,label = TeX("\\textbf{1.22% (0 - 30.5) virions fully infectious}")) +
ylab("% Fully Infectious Virions")
ggsave("Figures/2B_Percent_Fully_Infectious_Particles.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")
write.csv(Geom %>% mutate(Fully_Infectious = Particle_PFU * 100) %>% dplyr::select(Pp,Fully_Infectious),
file = file.path(Data.Path,"Fig2B_Data.csv"), row.names = FALSE)
# Figure 2C: p(Productive) vs. Virion Number ----
Fig.2 +
geom_line(data = Segment.Probs %>% filter(Segments == 8,
floor(Pp*100) %in% c(10,58,90)),
aes(x = Virion,
y = Prob * 100,
group = Pp,
color = Pp),
lwd = 1) +
labs(x="Virions Infecting",
y="Expected % Productive") +
scale_color_viridis(begin = 0.1, end = 0.9) +
scale_x_continuous(trans = "log10")
ggsave("Figures/2C_Percent_Cells_Productive_v_Virion_Number.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")
write.csv(Segment.Probs %>% filter(Segments == 8, floor(Pp*100) %in% c(10,58,90)) %>% dplyr::select(Pp,Virion,Prob),
file = file.path(Data.Path,"Fig2C_Data.csv"), row.names = FALSE)
# Figure 2D: Infectious Unit vs. Pp ----
Fig.2 +
geom_line(data = Geom,
aes(x = Pp,
y = Inf.Unit),
lwd = 1) +
labs(x=bquote(bold(P[P])),
y="Infectious Unit (# Virions)") +
scale_y_continuous(limits = c(0,55), breaks = c(0,10,20,30,40,50))## Warning: Removed 4 rows containing missing values (geom_path).
ggsave("Misc_Plots/Infectious_Unit_No_Data.pdf",
dpi = 300,
width = 4,
height = 4,
unit = "in")## Warning: Removed 4 rows containing missing values (geom_path).Fig.2 +
geom_line(data = Geom,
aes(x = Pp,
y = Inf.Unit),
lwd = 1) +
labs(x=bquote(bold(P[P])),
y="Infectious Unit (# Virions)") +
scale_y_continuous(breaks = c(0,3.6,10,20,30,40,50),
limits = c(0,55)) +
geom_ribbon(data = Geom %>% filter(Pp > 0.03),
aes(x = Pp,
ymax = pmin(Inf.Unit + Inf.Unit.sd * 1.96,55),
ymin = pmax(Inf.Unit - Inf.Unit.sd * 1.96,1)),
alpha = 0.3) +
geom_point(data = Inf.Unit.Sum,
aes(x = Pp,
y = Inf.Unit,
color = factor(Rep)),
size = 2) +
guides(color = FALSE) +
xlab(bquote(bold(P[P]))) +
scale_color_tableau(palette = "Tableau 20") +
geom_segment(aes(x = 0, xend = 0.52,
y = 3.6, yend = 3.6),
lty = 2) +
# annotate("text",x = 0.6, y = 40,size = 5,label = TeX("\\textbf{Infectious Unit = 3.6 (1 - 6.5) virions}")) +
ggsave("Figures/2D_Infectious_Unit_with_Data.pdf",
dpi = 300,
width = 5,
height = 5,
unit = "in")## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
write.csv(Geom %>% dplyr::select(Pp,Inf.Unit),
file = file.path(Data.Path,"Fig2D_Data.csv"), row.names = FALSE)
tau = matrix(nrow = 1,
ncol = 9,
data = c(1,0,0,0,0,0,0,0,0))
sub.tau = matrix(nrow = 1,
ncol = 8,
data = c(1,0,0,0,0,0,0,0))
trans = matrix(nrow = 9,
ncol = 9,
data = 0)
ident = diag(8)
one.mat = matrix(nrow = 8,
ncol = 1,
data = 1)
# Define Transition Matrix for Pp = 0.58
for (i in 0:8) {
for (j in i:8) {
trans[i+1,j+1] = dbinom(x = j - i, size = 8 - i, prob = 0.58)
}
trans[9,9] = 1
}
# Calculate Infectious Unit
sub.trans = trans[1:8,1:8]
Q = sub.trans
N = solve(ident - Q)
t = N %*% one.mat
tsq = t^2
Inf.Unit = sub.tau %*% solve(ident - sub.trans) %*% one.mat
Inf.Unit.sd = sqrt(((2 * N - ident) %*% t - tsq)[1,1])
Inf.Unit # 3.6## [,1]
## [1,] 3.633322 Inf.Unit + 1.96 * Inf.Unit.sd # 6.5## [,1]
## [1,] 6.480873 Inf.Unit - 1.96 * Inf.Unit.sd # 0.79## [,1]
## [1,] 0.785772# Import Data ----
Sim3=read.csv(file.path(Data.Path,"NTJ_Sim3_Summary_Data.csv"),sep=",",header=TRUE)
Sim3=subset(Sim3,MOI>0)
Sim3$Infected=as.numeric(Sim3$Infected)
# Output:
# % Productive vs. MOI -- how many is enough/what is the point of diminishing returns?
# Segments/Cell vs. MOI -- slower accumulation of segments at lower Pp
# % Populations infected vs. MOI -- What is the infectious dose (MOI) at the population level?
# Misc Plots ----
Sum.Sim3 = Sim3 %>%
mutate(MOI = (MOI)) %>%
group_by(MOI,Pp) %>%
summarise(Segments = mean(Segments),
Productive = mean(Productive),
Infected = mean(Infected) * 100,
Semi = mean(Semi),
Infected_Cells = mean(Infected_Cells),
Semi_Prop = mean(Semi) / mean(Infected_Cells))
Fig.3 = ggplot() +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5, color = "black"),
axis.line.y = element_line(size=0.5, color = "black"),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
Fig.3 +
geom_line(data = Sum.Sim3 %>% filter(Pp != 0.58),aes(x = log10(MOI),
y = Segments,
color = Pp,
group = Pp), lwd = 1) +
scale_color_viridis(guide = FALSE) +
labs(x = TeX("\\textbf{MOI (log_1_0 Virions / Cell)}"),
y = TeX("\\textbf{Segments / Cell}"),
color = NULL) +
ggsave('Misc_Plots/Population_Segments_Per_Cell.pdf',
width = 4,
height = 4,
unit = "in")
# % Cells Productive vs. MOI
Fig.3 +
geom_line(data=Sum.Sim3 %>% filter(Pp != 0.58),aes(x = log10(MOI),
y = Productive,
color = Pp,
group = Pp), lwd = 1) +
scale_color_viridis(guide = FALSE) +
labs(x = TeX("\\textbf{MOI (log_1_0 Virions / Cell)}"),
y = TeX("\\textbf{% Cells Productive}"),
color = NULL)
ggsave('Misc_Plots/Population_Percent_Cells_Productive.pdf',
width = 4,
height = 4,
unit = "in")
# Figure 3A: Populations Infected vs. MOI ----
Fig.3 +
geom_line(data=Sum.Sim3 %>% filter(Pp != 0.58),
aes(x = log10(MOI),
y = Infected,
color = Pp,
group = Pp),
lwd = 1) +
scale_color_viridis(guide = FALSE) +
labs(x = TeX("\\textbf{MOI (log_1_0 Virions / Cell)}"),
y = TeX("\\textbf{% Populations Infected}"),
color = NULL)
ggsave('Figures/3A_Populations_Infected.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(Sum.Sim3 %>% filter(Pp != 0.58) %>% dplyr::select(Pp,MOI,Infected),
file = file.path(Data.Path,"Fig3A_Data.csv"), row.names = FALSE)
# Figure 3B: Population ID50 vs. Pp ----
Pop.ID50 = matrix(nrow = 10,
ncol = 4,
data = 0)
colnames(Pop.ID50) = c("Pp","ID50","Low95","High95")
Pop.ID50[,"Pp"] = 1:10 / 10
for (p in 1:9) {
model = drm(data = Sum.Sim3 %>% filter(Pp == (p / 10), Infected < 100),
formula = Infected ~ MOI,
fct = LL.5()) %>%
ED(50)
Pop.ID50[p,"ID50"] = model[,1]
Pop.ID50[p,"Low95"] = model[,1] - 1.96 * model[,2]
Pop.ID50[p,"High95"] = model[,1] + 1.96 * model[,2]
}## Warning in sqrt(dEDval %*% varCov %*% dEDval): NaNs produced
## Warning in sqrt(dEDval %*% varCov %*% dEDval): NaNs produced##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 4.1287 NA
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.53831 0.34448
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.0439904 0.0005689
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.0088626 0.0001555
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 1.7516e-03 8.6606e-05
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 4.5202e-04 6.1038e-05
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 1.2799e-04 1.1764e-05
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 6.7929e-05 1.7795e-05
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 4.4033e-05 1.1480e-06Pop.ID50[10,"ID50"] = 1e-5
Pop.ID50 = data.frame(Pop.ID50)
Pop.ID50$ID50 = log10(Pop.ID50$ID50)
Pop.ID50$Low95 = log10(Pop.ID50$Low95)## Warning: NaNs producedPop.ID50$High95 = log10(Pop.ID50$High95)
Pop.ID50$Num_Full = 1 / (Pop.ID50$Pp ^ 8)
Pop.ID50$MOI_Full = (Pop.ID50$Num_Full / 1e5) %>% log10
Fig.3 +
geom_line(data = Pop.ID50[1:9,],
aes(x = Pp,
y = ID50),
lwd = 1) +
geom_line(data = Pop.ID50[1:9,],
aes(x = Pp,
y = MOI_Full),
lty = 2,
lwd = 1) +
geom_ribbon(data = Pop.ID50[1:9,],
aes(x = Pp,
ymin = Low95,
ymax = High95),
alpha = 0.3) +
scale_x_continuous(limits = c(0,1), breaks = c(0,0.25,0.5,0.75,1)) +
scale_y_continuous(limits = c(-5,3), breaks = c(-4,-2,0,2)) +
labs(x = TeX("\\textbf{P_P}"),
y = TeX("\\textbf{Population ID_5_0 (log_1_0 Virions / Cell)}")) +
theme(text = element_text(size=14,face="bold"))
ggsave('Figures/3B_Population_ID50.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(Pop.ID50 %>% rename(Complementation_ID50 = ID50, No_Complementation_ID50 = MOI_Full) %>% select(Pp, Complementation_ID50, No_Complementation_ID50),file = file.path(Data.Path,"Fig3B_Data.csv"), row.names = FALSE)
Fig.3 +
geom_line(data=Sum.Sim3,aes(x = log10(MOI),
y = Semi,
color = Pp,
group = Pp), lwd = 1.2) +
scale_color_viridis(guide = FALSE) +
labs(x = TeX("\\textbf{MOI (log_1_0 Virions / Cell)}"),
y = TeX("\\textbf{% Cells Semi-Infected}"),
color = NULL)
ggsave('Misc_Plots/Population_Percent_Cells_with_IVGs.pdf',
width = 4,
height = 4,
unit = "in")The following code chunk describes the process of summarizing the full simulation data from the individual-based model of viral replication in terms of peak Productive Cell, Virions, and Semi-Infected cell numbers, as well as the extraction of the data from each hour of the simulations, followed by the extraction of hour-by-hour data from each simulation for plotting the dynamics of infection in Supplementary Figure 2.
The full data set is 1.75 GB, and so was not uploaded to Github, but is available upon request. The script to generate the full data set can be found in the “Simulations” folder.
NTJ21 = read.csv(file = file.path(Data.Path,"NTJ_Sim21_Data.csv"))
NTJ21.Peak = NTJ21 %>%
group_by(Sim,Spread_Frac,Pp,Virus_D) %>%
summarise(Max_Productive = max(Productive),
Max_Virions = max(Virions),
Max_Semi = max(Infected - Productive))
NTJ21.Hours = NTJ21 %>%
filter((t * 10) %in% (1:(96 * 10)))
write.csv(NTJ21.Peak,file=file.path(Data.Path,"NTJ21_Peak_Data.csv"),row.names=FALSE)
write.csv(NTJ21.Hours,file=file.path(Data.Path,"NTJ21_Hours_Data.csv"),row.names=FALSE)# Import Data ----
NTJ21.Peak = read.csv(file = file.path(Data.Path,"NTJ21_Peak_Data.csv")) %>%
mutate(Virus_D = log10(Virus_D / 3600)) %>%
filter(Virus_D < 5,
Spread_Frac %in% c(0))
NTJ21.Hours = read.csv(file=file.path(Data.Path,"NTJ21_Hours_Data.csv")) %>%
mutate(Virus_D = log10(Virus_D / 3600)) %>%
filter(Virus_D < 5,
Spread_Frac %in% c(0))
NTJ21.Hours = NTJ21.Hours %>%
mutate(Spread_Frac = as.numeric(Spread_Frac))
NTJ21.Hours = NTJ21.Hours %>%
mutate(Spread_Frac = factor(Spread_Frac,levels = c("0.25","0")),
Spread_Frac = recode(Spread_Frac,"0.25" = "25% Local Spread","0" = "Diffusion Only"))
NTJ21.Peak = NTJ21.Peak %>%
mutate(Spread_Frac = factor(Spread_Frac,levels = c("0.25","0")),
Spread_Frac = recode(Spread_Frac,"0.25" = "25% Local Spread","0" = "Diffusion Only"))
# Plot Stock ----
Fig.4 = ggplot() +
scale_color_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1"), guide = FALSE) +
scale_fill_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1"), guide = FALSE) +
scale_y_continuous(limits = c(0,105)) +
geom_vline(xintercept = log10((5.825)), lty = 2) +
coord_cartesian(xlim = c(-0.5,4)) +
theme(text=element_text(size = 15,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
# Figure 4A: Initial Growth Rate ----
Early.Time = 12
df = NTJ21.Hours %>%
filter(t == Early.Time) %>%
mutate(Slope = pmax(log10(Productive),0)) %>%
select(Pp,Virus_D,Spread_Frac,Sim,Slope) %>%
group_by(Pp,Virus_D,Spread_Frac) %>%
summarise(r0 = mean(Slope)) %>%
ungroup %>%
mutate(Pp = Pp + 0.005) %>%
mutate(Pp = round(Pp,2))
Fig.4 +
geom_point(data = df,
aes(x = Virus_D,
y = r0,
color = Spread_Frac)) +
geom_smooth(data = df,
aes(x = Virus_D,
y = r0,
linetype = factor(Pp),
color = Spread_Frac,
fill = Spread_Frac,
group = interaction(Spread_Frac,Pp)),
lwd = 1.2,
alpha = 0.3) +
labs(y = TeX("\\textbf{$\\log_1_0$ Cells Infected in 12 Hours}"),
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}"),
linetype = TeX("\\textbf{$\\P_P$}"),
size = TeX("\\textbf{$\\P_P$}"),
color = NULL,
fill = NULL) +
scale_linetype_manual(values = c("1" = 1,"0.58" = 1212)) +
scale_y_continuous(limits = c(0,4)) +
scale_color_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1")) +
scale_fill_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1"))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave('Figures/4A_r0_Raw.pdf',
width = 6,
height = 6,
unit = "in")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'write.csv(df, file = file.path(Data.Path,"Fig4A_Data.csv"), row.names = FALSE)
write.csv(NTJ21.Hours %>%
filter(t == Early.Time) %>%
mutate(r0 = pmax(log10(Productive),0)) %>%
select(Pp,Virus_D,Spread_Frac,Sim,r0), file = file.path(Data.Path,"Fig4A_Full_Data.csv"), row.names = FALSE)
# Figure 4B: Reduction in Initial Growth Rate ----
Early.Time = 12
df1 = NTJ21.Hours %>%
filter(t == Early.Time,
Pp == 1) %>%
mutate(Slope = log10(Productive) / t) %>%
group_by(Virus_D,Spread_Frac) %>%
summarise(Intact_r0 = mean(Slope))
df2 = NTJ21.Hours %>%
filter(t == Early.Time,
Pp < 1) %>%
mutate(Slope = pmax(log10(Productive) / t,0)) %>%
select(Pp,Virus_D,Spread_Frac,Sim,Slope)
df3 = right_join(df1,df2) %>%
mutate(Cost = (1 - (Slope / Intact_r0)) * 100) %>%
group_by(Virus_D,Spread_Frac) %>%
summarise(Cost = mean(Cost))## Joining, by = c("Virus_D", "Spread_Frac")Fig.4 +
geom_point(data = df3,
aes(x = Virus_D,
y = Cost,
color = Spread_Frac,
group = Spread_Frac)) +
geom_smooth(data = df3,
aes(x = Virus_D,
y = Cost,
group = Spread_Frac,
color = Spread_Frac,
fill = Spread_Frac),
lty = 1212,
alpha = 0.3) +
labs(y = "% Reduction in Initial Growth Rate",
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}"))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave('Figures/4B_r0_Reduction.pdf',
width = 6,
height = 6,
unit = "in")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'write.csv(df3, file = file.path(Data.Path,"Fig4B_Data.csv"), row.names = FALSE)
write.csv(right_join(df1,df2) %>%
mutate(Cost = (1 - (Slope / Intact_r0)) * 100), file = file.path(Data.Path,"Fig4B_Full_Data.csv"), row.names = FALSE)## Joining, by = c("Virus_D", "Spread_Frac")# Figure 4C: Increase in Time to Infect 100 Cells ----
Target_Cells = 100
df1 = NTJ21.Hours %>%
filter(Pp == 1,
Productive >= Target_Cells) %>%
group_by(Spread_Frac,Virus_D,Sim) %>%
summarise(To_Target_Cells = min(t)) %>%
ungroup() %>%
group_by(Spread_Frac,Virus_D) %>%
summarise(Intact_To_Target_Cells = mean(To_Target_Cells))
df2 = NTJ21.Hours %>%
filter(Pp < 1,
Productive >= Target_Cells) %>%
group_by(Spread_Frac,Virus_D,Sim) %>%
summarise(To_Target_Cells = min(t)) %>%
select(Spread_Frac,Virus_D,Sim,To_Target_Cells)
df3 = right_join(df1,df2) %>%
mutate(Target_Cell_Cost = (To_Target_Cells - Intact_To_Target_Cells) / Intact_To_Target_Cells * 100) %>%
group_by(Virus_D,Spread_Frac) %>%
summarise(Target_Cell_Cost = mean(Target_Cell_Cost))## Joining, by = c("Spread_Frac", "Virus_D")Fig.4 +
geom_point(data = df3,
aes(x = Virus_D,
y = Target_Cell_Cost,
color = Spread_Frac,
group = Spread_Frac)) +
geom_smooth(data = df3,
aes(x = Virus_D,
y = Target_Cell_Cost,
group = Spread_Frac,
color = Spread_Frac,
fill = Spread_Frac),
lty = 1212,
alpha = 0.3) +
scale_y_continuous(limits = c(0,600)) +
labs(y = "% Increase in Time to Infect 100 Cells",
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}"))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave('Figures/4C_Increase_in_Time_to_100_Cells.pdf',
width = 6,
height = 6,
unit = "in")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'write.csv(df3, file = file.path(Data.Path,"Fig4C_Data.csv"), row.names = FALSE)
write.csv(right_join(df1,df2) %>%
mutate(Target_Cell_Cost = (To_Target_Cells - Intact_To_Target_Cells) / Intact_To_Target_Cells * 100), file = file.path(Data.Path,"Fig4C_Full_Data.csv"), row.names = FALSE)## Joining, by = c("Spread_Frac", "Virus_D")# Figure 4D: Increase in Time to Produce 10^5 Virions ----
Target_Virions = 1e5
df1 = NTJ21.Hours %>%
filter(Pp == 1,
Virions >= Target_Virions) %>%
group_by(Spread_Frac,Virus_D,Sim) %>%
summarise(To_Target_Virions = min(t)) %>%
ungroup() %>%
group_by(Spread_Frac,Virus_D) %>%
summarise(Intact_To_Target_Virions = mean(To_Target_Virions))
df2 = NTJ21.Hours %>%
filter(Pp < 1,
Virions >= Target_Virions) %>%
group_by(Spread_Frac,Virus_D,Sim) %>%
summarise(To_Target_Virions = min(t)) %>%
select(Spread_Frac,Virus_D,Sim,To_Target_Virions)
df3 = right_join(df1,df2) %>%
mutate(Target_Virion_Cost = (To_Target_Virions - Intact_To_Target_Virions) / Intact_To_Target_Virions * 100) %>%
group_by(Virus_D,Spread_Frac) %>%
summarise(Target_Virion_Cost = mean(Target_Virion_Cost))## Joining, by = c("Spread_Frac", "Virus_D")Fig.4 +
geom_point(data = df3,
aes(x = Virus_D,
y = Target_Virion_Cost,
color = Spread_Frac,
group = Spread_Frac)) +
geom_smooth(data = df3,
aes(x = Virus_D,
y = Target_Virion_Cost,
group = Spread_Frac,
color = Spread_Frac,
fill = Spread_Frac),
lty = 1212,
alpha = 0.3) +
scale_y_continuous(limits = c(0,400)) +
labs(y = TeX("\\textbf{% Increase in Time to Yield 10^5 Virions}"),
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}"))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave('Figures/4D_Increase_in_Time_to_1e5_Virions.pdf',
width = 6,
height = 6,
unit = "in")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'write.csv(df3, file = file.path(Data.Path,"Fig4D_Data.csv"), row.names = FALSE)
write.csv(right_join(df1,df2) %>%
mutate(Target_Virion_Cost = (To_Target_Virions - Intact_To_Target_Virions) / Intact_To_Target_Virions * 100), file = file.path(Data.Path,"Fig4D_Full_Data.csv"), row.names = FALSE)## Joining, by = c("Spread_Frac", "Virus_D")# Supplemental Figure 2: Representative Time Course of Productive Cells vs. Time ----
ggplot() +
geom_line(data = NTJ21.Hours %>% filter(Pp == 1,Sim == 1,Spread_Frac == "Diffusion Only"),
aes(x = t,
y = Productive,
color = Virus_D,
group = Virus_D),
lwd = 1.2) +
scale_color_viridis(guide = FALSE) +
scale_y_continuous(limits = c(0,10000),labels = comma) +
scale_x_continuous(breaks = c(0,24,48,72)) +
labs(x = TeX("\\textbf{Time Post-Inoculation (Hours)}"),
y = TeX("\\textbf{Productively Infected Cells}")) +
theme(text=element_text(size = 13,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave('Figures/Supp2A_Time_Course_Pp1.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ21.Hours %>% filter(Pp == 1,Sim == 1,Spread_Frac == "Diffusion Only") %>% select(t,Pp,Virus_D,Productive), file = file.path(Data.Path,"Supp2A_Data.csv"), row.names = FALSE)
ggplot() +
geom_line(data = NTJ21.Hours %>% filter(Pp == 0.575,Sim == 1,Spread_Frac == "Diffusion Only"),
aes(x = t,
y = Productive,
color = Virus_D,
group = Virus_D),
lwd = 1.2) +
scale_color_viridis(guide = FALSE) +
scale_y_continuous(limits = c(0,10000),labels = comma) +
scale_x_continuous(breaks = c(0,24,48,72)) +
labs(x = TeX("\\textbf{Time Post-Inoculation (Hours)}"),
y = TeX("\\textbf{Productively Infected Cells}")) +
theme(text=element_text(size = 13,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave('Figures/Supp2C_Time_Course_Pp_058_.pdf',
width = 6,
height = 5,
unit = "in")
write.csv(NTJ21.Hours %>% filter(Pp == 0.575,Sim == 1,Spread_Frac == "Diffusion Only") %>% select(t,Pp,Virus_D,Productive), file = file.path(Data.Path,"Supp2C_Full_Data.csv"), row.names = FALSE)
# Supplemental Figure 2: Representative Time Courses of Virions vs. Time ----
ggplot() +
geom_line(data = NTJ21.Hours %>% filter(Pp == 1,Sim == 1,Spread_Frac == "Diffusion Only"),
aes(x = t,
y = Virions %>% log10,
color = Virus_D,
group = Virus_D),
lwd = 1.2) +
scale_color_viridis(guide = FALSE) +
scale_x_continuous(breaks = c(0,24,48,72)) +
labs(x = TeX("\\textbf{Time Post-Inoculation (Hours)}"),
y = TeX("\\textbf{$\\log_1_0$ Virions}")) +
theme(text=element_text(size = 13,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave('Figures/Supp2B_Virion_Time_Course_PP1.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ21.Hours %>% filter(Pp == 1,Sim == 1,Spread_Frac == "Diffusion Only") %>% select(t,Pp,Virus_D,Virions), file = file.path(Data.Path,"Supp2B_Data.csv"), row.names = FALSE)
ggplot() +
geom_line(data = NTJ21.Hours %>% filter(Pp == 0.575,Sim == 1,Spread_Frac == "Diffusion Only"),
aes(x = t,
y = Virions %>% log10,
color = Virus_D,
group = Virus_D),
lwd = 1.2) +
scale_color_viridis(guide = FALSE) +
scale_x_continuous(breaks = c(0,24,48,72)) +
labs(x = TeX("\\textbf{Time Post-Inoculation (Hours)}"),
y = TeX("\\textbf{$\\log_1_0$ Virions}")) +
theme(text=element_text(size = 13,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave('Figures/Supp2D_Virion_Time_Course_Pp_058.pdf',
width = 6,
height = 5,
unit = "in")
write.csv(NTJ21.Hours %>% filter(Pp == 0.575,Sim == 1,Spread_Frac == "Diffusion Only") %>% select(t,Pp,Virus_D,Virions), file = file.path(Data.Path,"Supp2D_Data.csv"), row.names = FALSE)
D_Legend = ggplot() +
geom_line(data = NTJ21.Hours %>% filter(Pp == 0.575,Sim == 1,Spread_Frac == "Diffusion Only"),
aes(x = t,
y = Virions %>% log10,
color = Virus_D,
group = Virus_D),
lwd = 1.2) +
scale_color_viridis() +
scale_x_continuous(breaks = c(0,24,48,72)) +
labs(color = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}")) +
theme(text=element_text(size = 13,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
# Supplemental Figure 2 Legend: Diffusion Coefficient Color Gradient ----
plot(g_legend(D_Legend))
ggsave('Figures/Supp2_D_Legend.pdf',
width = 4,
height = 4,
unit = "in")
Supp.Fig.2.Legend = arrangeGrob(g_legend(D_Legend))
ggsave(Supp.Fig.2.Legend,
file = "Figures/Supp2Legend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")
Legend.Fig.4 = ggplot() +
geom_point(data = df,
aes(x = Virus_D,
y = r0,
color = Spread_Frac)) +
geom_smooth(data = df,
aes(x = Virus_D,
y = r0,
linetype = factor(Pp),
color = Spread_Frac,
fill = Spread_Frac,
group = interaction(Spread_Frac,Pp)),
lwd = 1.2,
alpha = 0.3) +
labs(y = TeX("\\textbf{$\\log_1_0$ Cells Infected in 12 Hours}"),
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}"),
linetype = TeX("\\textbf{$\\P_P$}"),
size = TeX("\\textbf{$\\P_P$}"),
color = NULL,
fill = NULL) +
scale_linetype_manual(values = c("1" = 1,"0.58" = 1212)) +
scale_y_continuous(limits = c(0,4)) +
scale_color_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1")) +
scale_fill_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1")) +
theme(text=element_text(size = 15,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle = 0,vjust=0.75,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5),
axis.ticks.y = element_line(size = 0.5),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
Fig.4.Legend = arrangeGrob(g_legend(Legend.Fig.4))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'ggsave(Fig.4.Legend,
file = "Figures/4Legend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")# Import Data ----
Cell.Num = 6.4e5
NTJ29 = read.csv(file = file.path(Data.Path,"NTJ29_Data.csv"), header = TRUE) %>%
melt(id = c("MOI","Time")) %>%
na.omit() %>%
mutate(PFU = (value))
# Analyze Data ----
# Growth Curve Summary
NTJ29.Sum = NTJ29 %>% mutate(PFU = log10(PFU),
MOI = (MOI)) %>%
group_by(MOI, Time) %>%
summarise(Titer = mean(PFU),
se = sd(PFU)/sqrt(length(PFU)))
# Virus Production by Time Window
NTJ29.Kinetics = NTJ29 %>% group_by(MOI) %>%
mutate(growth = 2 * (PFU - lag(PFU))) %>%
mutate(growth = pmax(growth,1)) %>%
filter(Time %in% c(12,16,24)) %>%
mutate(growth = log10(growth)) %>%
group_by(MOI,Time) %>%
summarise(Delta = mean(growth),
se = sd(growth)/sqrt(length(growth)))
# HA expression by flow cytometry (12 hours)
NTJ29.Flow = read.csv(file = file.path(Data.Path,"NTJ29_Flow_Data.csv"), header = TRUE) %>%
group_by(MOI) %>%
summarise(HA.Prop = mean(HA_Cells),
HA.Prop.se = sd(HA_Cells)/sqrt(length(HA_Cells)),
HA.MFI = mean(HA_MFI),
HA.MFI.se = sd(HA_MFI)/sqrt(length(HA_MFI)))
# Burst Size
NTJ29.Burst = NTJ29 %>%
filter(Time >= 12) %>%
mutate(Total.PFU = PFU * 2,
Burst = Total.PFU / ((1 - ppois(0, lambda = MOI, lower = TRUE)) * Cell.Num),
Amp = Total.PFU / (Cell.Num * MOI)) %>%
group_by(MOI,Time) %>%
summarise(se.Burst = sd(Burst)/sqrt(length(Burst)),
Burst = mean(Burst),
se.Amp = sd(Amp)/sqrt(length(Amp)),
Amp = mean(Amp))
# Stock Plot ----
Fig.5 = ggplot() +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5, color = "black"),
axis.line.y = element_line(size=0.5, color = "black"),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
# Figure 5A: Growth Curve ----
NTJ29.Growth.Plot = Fig.5 +
geom_line(data = NTJ29.Sum %>% filter(MOI >= 1),
aes(x = Time,
y = Titer,
group = factor(MOI),
color = factor(MOI)),
lwd = 1) +
geom_errorbar(data = NTJ29.Sum %>% filter(MOI >= 1), aes(x = Time, group = MOI, ymin = Titer - se, ymax = Titer + se),
color = "black",
width = 1,
lwd = .75) +
labs(x = "Time Post-Inoculation (Hours)",
y = TeX("\\textbf{Viral Titers ($\\log_1_0$ PFU/mL)}"),
color = "MOI (PFU/cell)") +
scale_color_viridis(discrete = TRUE) +
scale_x_continuous(limits = c(0,50), breaks = c(0,12,24,48)) +
geom_hline(yintercept = log10(50),lty = 2)
# NTJ29.Legend = arrangeGrob(g_legend(NTJ29.Growth.Plot +theme(text = element_text(size = 12, face = "bold"))))
#
# ggsave(NTJ29.Legend,
# file = "Figures/5Legend.pdf",
# dpi = 300,
# width = 5,
# height = 4,
# unit = "in")
NTJ29.Growth.Plot + theme(legend.position = "NA")
ggsave('Figures/5A_PFU_Growth_Curve.pdf',
dpi = 300,
width = 5,
height = 4,
unit = "in")
write.csv(NTJ29.Sum, file = file.path(Data.Path,"Fig5A_Data.csv"), row.names = FALSE)
# Figure 5B: %HA+ Cells vs. MOI ----
Fig.5 +
geom_line(data = NTJ29.Flow %>% filter(MOI >= 0),
aes(x = MOI,
y = HA.Prop),
lwd = 1) +
geom_errorbar(data = NTJ29.Flow %>% filter(MOI >= 0), aes(x = MOI, group = MOI, ymin = HA.Prop - HA.Prop.se, ymax = HA.Prop + HA.Prop.se),
color = "black",
width = .5,
lwd = .75) +
labs(x = "Multiplicity of Infection (PFU / cell)",
y = TeX("\\textbf{% Cells Expressing HA^+}")) +
scale_y_continuous(limits = c(0,100), breaks = c(0,20,40,60,80,100)) +
scale_x_continuous(limits = c(0,20), breaks = c(1,3,6,10,20))
ggsave('Figures/5B_Percent_HA_Positive.pdf',
dpi = 300,
width = 5,
height = 4,
unit = "in")
write.csv(NTJ29.Flow, file = file.path(Data.Path,"Fig5B_Data.csv"), row.names = FALSE)
# Figure 5C: Burst Size vs. MOI ----
Fig.5 +
geom_bar(data = NTJ29.Burst %>% filter(Time == 48, MOI >= 1),
stat = "identity",
aes(x = factor(MOI),
y = Burst,
group = MOI,
fill = (MOI)),
lwd = 1) +
scale_fill_viridis(trans = "log10") +
# annotate("text",x = 2.5, y = 14,size = 5,label = TeX("\\textbf{Burst Size = 11.5 (10.6 - 12.5) PFU/cell}")) +
geom_errorbar(data = NTJ29.Burst %>% filter(Time == 48, MOI >= 1), aes(x = factor(MOI), group = MOI, ymin = Burst - se.Burst, ymax = Burst + se.Burst),
color = "black",
width = .3,
lwd = .75) +
labs(x = "Multiplicity of Infection (PFU / cell)",
y = TeX("\\textbf{Burst Size (PFU / cell)}"),
fill = "MOI") +
scale_color_viridis(trans = "log10") +
scale_y_continuous(limits = c(0,15), breaks = c(0,5,10,15,20))
ggsave('Figures/5C_Burst_Size.pdf',
dpi = 300,
width = 5,
height = 4,
unit = "in")
write.csv(NTJ29.Burst %>% filter(Time == 48, MOI >= 1), file = file.path(Data.Path,"Fig5C_Data.csv"), row.names = FALSE)
# Supp Figure 3A: Virus production by time window ----
ggplot() +
geom_bar(data = NTJ29.Kinetics %>% filter(Time <= 24),
stat = "identity",
position = "dodge",
aes(x = factor(MOI),
y = Delta,
fill = MOI)) +
geom_errorbar(data = NTJ29.Kinetics,
position = "dodge",
aes(x = factor(MOI), ymin = Delta - se, ymax = Delta + se),
color = "black",
width = .3,
lwd = .75) +
facet_grid(.~Time,labeller = labeller(Time = c("12" = "8 h -> 12 h",
"16" = "12 h -> 16 h",
"24" = "16 h -> 24 h"))) +
scale_fill_viridis(trans = "log10") +
labs(x = "Multiplicity of Infection (PFU / cell)",
y = TeX("\\textbf{Virus Production ($\\log_1_0$ PFU)}")) +
theme(text=element_text(size=16,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks = element_line(size=0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
ggsave('Figures/Supp3A_Delta_PFU_Bars.pdf',
dpi = 300,
width = 8,
height = 6,
unit = "in")
write.csv(NTJ29.Kinetics %>% filter(Time <= 24), file = file.path(Data.Path,"Supp3A_Data.csv"), row.names = FALSE)
# Supp Figure 3B: HA Intensity vs. MOI ----
ggplot() +
geom_line(data = NTJ29.Flow,
aes(x = MOI,
y = HA.MFI),
lwd = 1) +
geom_errorbar(data = NTJ29.Flow, aes(x = MOI, group = MOI, ymin = HA.MFI - HA.MFI.se, ymax = HA.MFI + HA.MFI.se),
color = "black",
width = .5,
lwd = .75) +
labs(x = "Multiplicity of Infection (PFU / cell)",
y = TeX("\\textbf{HA Expression (Median FI)}")) +
scale_y_continuous(limits = c(0,20000), breaks = c(1000,5000,10000,15000,20000), labels = comma) +
scale_x_continuous(limits = c(0,21), breaks = c(1,3,6,10,20)) +
theme(text=element_text(size=12,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks = element_line(size=0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave('Figures/Supp3B_HA_MFI.pdf',
dpi = 300,
width = 4,
height = 3,
unit = "in")
write.csv(NTJ29.Flow, file = file.path(Data.Path,"Supp3B_Data.csv"), row.names = FALSE)
# Supp Figure 3C: Virus Amplification vs. MOI ----
ggplot() +
geom_bar(data = NTJ29.Burst %>% filter(Time == 48, MOI >= 1),
stat = "identity",
aes(x = factor(MOI),
y = Amp,
group = MOI,
fill = MOI),
lwd = 1) +
scale_fill_viridis(trans = "log10") +
geom_errorbar(data = NTJ29.Burst %>% filter(Time == 48, MOI >= 1), aes(x = factor(MOI), group = MOI, ymin = Amp - se.Amp, ymax = Amp + se.Amp),
color = "black",
width = .3,
lwd = .75) +
geom_hline(yintercept = 1,lty = 2) +
labs(x = "Multiplicity of Infection (PFU / cell)",
y = TeX("\\textbf{Virus Amplification (PFU output / input)}"),
fill = "MOI") +
scale_color_viridis(trans = "log10") +
scale_y_continuous(limits = c(0,7), breaks = c(0,1,2,3,4,5,6,7,8,10)) +
theme(text=element_text(size=12,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks = element_line(size=0.5, color = "black"),
axis.title.y = element_text(size=rel(1.0),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
ggsave('Figures/Supp3C_Virus_Amplification.pdf',
dpi = 300,
width = 4.5,
height = 3.6,
unit = "in")
write.csv(NTJ29.Burst %>% filter(Time == 48, MOI >= 1), file = file.path(Data.Path,"Supp3C_Data.csv"), row.names = FALSE)
# Calculate CI for burst size -----
NTJ29.Burst.Raw = NTJ29 %>%
filter(Time >= 12) %>%
mutate(Total.PFU = PFU * 2,
Burst = Total.PFU / ((1 - ppois(0, lambda = MOI, lower = TRUE)) * Cell.Num),
Amp = Total.PFU / (Cell.Num * MOI))
NTJ29.Burst.Raw %>% filter(Time == 48,
MOI > 1) %>%
summarise(Burst.Mean = mean(Burst),
Burst.se = sd(Burst)/sqrt(length(Burst))) %>%
mutate(Burst.Max = Burst.Mean + 1.96 * Burst.se,
Burst.Min = Burst.Mean - 1.96 * Burst.se)## Burst.Mean Burst.se Burst.Max Burst.Min
## 1 11.5316 0.5003967 12.51237 10.55082# 11.5 (10.6 - 12.5)Fig.6 = ggplot() +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
# Figure 6A: Percent Infected Cells Containing IVGs vs. MOI ----
Fig.6 +
geom_line(data = Sum.Sim3 %>% filter(Pp %in% c(0.58)),
aes(x = log10(MOI),
y = Semi_Prop * 100,
group = Pp), lwd = 1.2, color = "black") +
scale_color_viridis(guide = FALSE, begin = 0.1, end = 0.9) +
labs(x = TeX("\\textbf{MOI (log_1_0 Virions / Cell)}"),
y = TeX("\\textbf{% Infected Cells with IVGs}")) +
scale_x_continuous(limits = c(-2,2)) ## Warning: Removed 22 rows containing missing values (geom_path).
ggsave('Figures/6A_Proportion_Infected_Cells_with_IVGs.pdf',
width = 5,
height = 5,
unit = "in")## Warning: Removed 22 rows containing missing values (geom_path).write.csv(Sum.Sim3 %>% filter(Pp %in% c(0.58)),
file = file.path(Data.Path,"Fig6A_Data.csv"), row.names = FALSE)
# Figure 6B: Peak # Cells with IVGs vs. Diffusion Coefficient ----
df1 = NTJ21.Peak %>%
group_by(Virus_D,Spread_Frac) %>%
summarise(Max_Semi = mean(Max_Semi))
Fig.6 +
geom_point(data = df1,
aes(x = Virus_D,
y = Max_Semi,
color = Spread_Frac)) +
geom_smooth(data = df1,
aes(x = Virus_D,
y = Max_Semi,
group = Spread_Frac,
color = Spread_Frac,
fill = Spread_Frac),
lty = 1212,
alpha = 0.3) +
scale_y_continuous(limits = c(-100,3200), breaks = c(0,1000,2000,3000), labels = comma) +
labs(y = TeX("\\textbf{Peak # Cells with IVGs}"),
x = TeX("\\textbf{Diffusion Coefficient ($\\log_1_0$ $\\mu$m^2/s)}")) +
scale_color_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1"), guide = FALSE) +
scale_fill_manual(values = c("Diffusion Only" = "midnightblue","25% Local Spread" = "steelblue1"), guide = FALSE) +
geom_vline(xintercept = log10((5.825)), lty = 2) +
coord_cartesian(xlim = c(-0.5,4)) +
theme(legend.position = "NA")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggsave('Figures/6B_Max_Semi.pdf',
width = 5,
height = 5,
unit = "in")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'write.csv(df1,
file = file.path(Data.Path,"Fig6B_Data.csv"), row.names = FALSE)
# Figure 6D: Enrichment vs. % WT HA+ Cells ----
NTJ52 = read.csv(file = file.path(Data.Path,"NTJ52_Data.csv"), header = TRUE) %>%
mutate(Expt = "NTJ52")
NTJ53 = read.csv(file = file.path(Data.Path,"NTJ53_Data.csv"), header = TRUE) %>%
mutate(Expt = "NTJ53")
NTJ54 = read.csv(file = file.path(Data.Path,"NTJ54_Data.csv"), header = TRUE) %>%
mutate(Expt = "NTJ54")
NTJ54.Meta = rbind(NTJ52,NTJ53,NTJ54) %>%
mutate(var_HA = Both + HA,
WT_HA = Both + His,
inv_WT_HA = 1 / WT_HA,
Both_HA = Both,
WT_Alone = His / (His + Neg),
WT_Co = Both / (Both + HA),
Enrichment = ((WT_Co) - (WT_Alone)) / (WT_Alone) * 100)
m1 = lm(data = NTJ54.Meta %>% filter(var_MOI * WT_MOI > 0),
formula = Enrichment ~ inv_WT_HA * Spread)
summary(m1)##
## Call:
## lm(formula = Enrichment ~ inv_WT_HA * Spread, data = NTJ54.Meta %>%
## filter(var_MOI * WT_MOI > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.545 -20.676 2.053 14.141 109.120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -343.24 22.20 -15.461 < 2e-16 ***
## inv_WT_HA 31327.60 1439.47 21.763 < 2e-16 ***
## Spread 294.62 31.07 9.482 9.34e-13 ***
## inv_WT_HA:Spread -25701.71 1859.28 -13.823 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.88 on 50 degrees of freedom
## Multiple R-squared: 0.9189, Adjusted R-squared: 0.914
## F-statistic: 188.8 on 3 and 50 DF, p-value: < 2.2e-16newdata = matrix(ncol = 4,
nrow = length(40:90)*2)
colnames(newdata) = c("WT_HA","inv_WT_HA","Spread","Enrichment")
newdata[,"WT_HA"] = rep(seq(from = 40, to = 90, by = 1),2)
newdata[,"inv_WT_HA"] = 1/rep(seq(from = 40, to = 90, by = 1),2)
newdata[,"Spread"] = c(rep(0,length(40:90)),rep(1,length(40:90)))
newdata = as.data.frame(newdata)
newdata$Enrichment = predict(m1,newdata,type="response",se.fit=TRUE)$fit
newdata$se = predict(m1,newdata,type="response",se.fit=TRUE)$se.fit
Fig.6 +
geom_point(data = NTJ54.Meta %>% filter(WT_MOI * var_MOI > 0),
aes(x = WT_HA,
y = Enrichment,
fill = factor(Spread),
color = factor(Spread)),
shape = 21,
size = 4,
stroke = 2) +
geom_line(data = newdata,
aes(x = WT_HA,
y = Enrichment,
group = Spread,
color = factor(Spread)),
lwd = 1.5) +
geom_ribbon(data = newdata %>% filter(Spread == 0),
aes(x = WT_HA,
y = Enrichment,
ymin = Enrichment - 1.96 * se,
ymax = Enrichment + 1.96 * se,
group = Spread),
fill = "navy",
alpha = 0.3) +
geom_ribbon(data = newdata %>% filter(Spread == 1),
aes(x = WT_HA,
y = Enrichment,
ymin = Enrichment - 1.96 * se,
ymax = Enrichment + 1.96 * se,
group = Spread),
fill = "green3",
alpha = 0.3) +
scale_x_continuous(limits = c(30,100)) +
scale_y_continuous(limits = c(-50, 450),breaks = c(0,100,200,300,400)) +
coord_cartesian(ylim = c(0,450)) +
scale_color_manual(values = c("0" = "navy","1" = "green3"),
guide = FALSE) +
scale_fill_manual(values = c("0" = NA, "1" = "green3"),
guide = FALSE) +
labs(x = TeX("\\textbf{% Cells WT HA^+}"),
y = "% WT Enrichment by Helper Virus")## Warning: Ignoring unknown aesthetics: y## Warning: Ignoring unknown aesthetics: y
ggsave('Figures/6D_Spatial_Structure_and_Enrichment.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ54.Meta %>% filter(WT_MOI * var_MOI > 0),
file = file.path(Data.Path,"Fig6D_Data.csv"), row.names = FALSE)
# Supplemental Figure 4C: Fold Change in % WT HA + vs. var MOI ----
NTJ51 = read.csv(file = file.path(Data.Path,"NTJ51_Data.csv"), header = TRUE) %>%
mutate(var_HA = Both + HA,
WT_HA = Both + His,
Both_HA = Both,
WT_Alone = His / (His + Neg),
WT_Co = Both / (Both + HA),
Enrichment = ((WT_Co) - (WT_Alone)) / (WT_Alone) * 100)
WT_base = NTJ51 %>% filter(MOI == 0) %>% dplyr::select(WT_HA) %>% sum() / 3
NTJ51 = NTJ51 %>% mutate(WT_add = WT_HA - WT_base,
WT_amp = WT_HA / WT_base)
Fig.6 +
geom_jitter(data = NTJ51 %>% filter(MOI > 0),
aes(x = MOI,
y = WT_amp)) +
geom_smooth(data = NTJ51,
aes(x = MOI,
y = WT_amp)) +
scale_x_continuous(trans = "log10") +
labs(y = "Fold Change %HA+ Cells (vs. No Help)",
x = "Helper MOI (PFU / cell)") +
geom_hline(yintercept = 1, lty = 2, lwd = 1.2) +
ggsave('Figures/Supp4C_WT_HA_Fold_Change_from_Helper.pdf',
width = 5,
height = 5,
unit = "in")## Warning: Transformation introduced infinite values in continuous x-axis## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Transformation introduced infinite values in continuous x-axis## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).
write.csv(NTJ51 %>% filter(MOI > 0),
file = file.path(Data.Path,"Supp4C_Data.csv"), row.names = FALSE)
# Supplemental Figure 4D: % VAR HA+ vs. % WT HA+ ----
Fig.6 +
geom_point(data = NTJ54.Meta %>% filter(var_MOI > 0),
size = 3,
stroke = 2,
aes(x = WT_HA,
y = var_HA,
color = factor(Spread),
shape = Expt,
fill = factor(Spread))) +
labs(x = TeX("\\textbf{% WT HA^+}"),
y = TeX("\\textbf{% Helper HA^+}"),
fill = "2° Spread",
color = "2° Spread") +
scale_shape_manual(values = c("NTJ52" = 21,
"NTJ53" = 22,
"NTJ54" = 24),
labels = c("NTJ52" = "1",
"NTJ53" = "2",
"NTJ54" = "3")) +
scale_color_manual(values = c("0" = "midnightblue", "1" = "green3")) +
scale_fill_manual(values = c("0" = NA, "1" = "green3"),guide = FALSE) +
scale_x_continuous(limits = c(0,100))
ggsave('Figures/Supp4D_Helper_v_WT_HA_in_Co_Infected_Cells.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ54.Meta %>% filter(var_MOI > 0),
file = file.path(Data.Path,"Supp4D_Data.csv"), row.names = FALSE)
var_M = lmer(data = NTJ54.Meta %>% filter(var_MOI > 0,
WT_MOI > 0),
formula = var_HA ~ WT_HA + Spread + (1|Expt))
summary(var_M)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: var_HA ~ WT_HA + Spread + (1 | Expt)
## Data: NTJ54.Meta %>% filter(var_MOI > 0, WT_MOI > 0)
##
## REML criterion at convergence: 394.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9613 -0.6918 -0.1281 0.7972 2.1330
##
## Random effects:
## Groups Name Variance Std.Dev.
## Expt (Intercept) 105.58 10.28
## Residual 86.11 9.28
## Number of obs: 54, groups: Expt, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 80.61456 7.95189 5.73291 10.138 7.11e-05 ***
## WT_HA -0.41479 0.07192 49.27879 -5.767 5.24e-07 ***
## Spread 10.34003 2.69630 48.99844 3.835 0.000359 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WT_HA
## WT_HA -0.637
## Spread -0.184 0.114t.test(NTJ54.Meta %>% filter(var_MOI > 0, WT_MOI > 0, Spread == 1) %>% dplyr::select(var_HA),
NTJ54.Meta %>% filter(var_MOI > 0, WT_MOI > 0, Spread == 0) %>% dplyr::select(var_HA))##
## Welch Two Sample t-test
##
## data: NTJ54.Meta %>% filter(var_MOI > 0, WT_MOI > 0, Spread == 1) %>% and NTJ54.Meta %>% filter(var_MOI > 0, WT_MOI > 0, Spread == 0) %>% dplyr::select(var_HA) and dplyr::select(var_HA)
## t = 3.3156, df = 41.204, p-value = 0.001915
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.734969 19.485587
## sample estimates:
## mean of x mean of y
## 63.51778 51.40750# Figure 7B: Measurement of RNA copies / mL M.STOP Virus Stock ----
NTJ60 = read.csv(file = file.path(Data.Path,"NTJ60_Data.csv"), header = TRUE)
NTJ60.ddPCR = NTJ60 %>% filter(Assay == "ddPCR",
Type %in% c("Mwt","M2stop","NS")) %>%
dcast(Strain + Rep ~ Type) %>%
mutate(Mtotal = M2stop + Mwt) %>%
melt(id = c("Strain","Rep")) %>%
na.omit %>%
mutate(Locus = str_replace_all(variable,c("Mwt" = "M", "M2stop" = "M")),
variable = str_replace_all(variable,c("Mwt" = "M2.Only", "M2stop" = "M1.Only")))## Using Readout as value column: use value.var to override.ggplot() +
geom_col(data = NTJ60.ddPCR %>%
filter(variable != "Mtotal",
Strain == "MUT") %>%
mutate(value = value * 20 / 4 * 40 * 20 / 12 * 60 / 280 * 1000 / 1e7),
aes(x = Locus,
y = value,
fill = variable),
color = "black",
position = "stack") +
facet_wrap(~Rep) +
scale_y_continuous(limits = c(0,10),breaks = 2 * 1:5, labels = comma) +
labs(x = NULL,
y = TeX("\\textbf{vRNA Concentration (copies * $10^7$ / mL)}"),
fill = "Segment") +
scale_fill_manual(values = c("M1.Only" = "steelblue",
"M2.Only" = "lightgoldenrod",
"NS" = "gray40")) +
theme(text=element_text(size=14,face="bold"),
strip.text.x = element_blank(),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),face = "bold",color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
ggsave('Figures/7B_Segment_Copy_in_Pan99_M_STOP_Stocks.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ60.ddPCR %>%
filter(variable != "Mtotal",
Strain == "MUT") %>%
mutate(value = value * 20 / 4 * 40 * 20 / 12 * 60 / 280 * 1000 / 1e7),
file = file.path(Data.Path,"Fig7B_Data.csv"), row.names = FALSE)
Legend.7B = ggplot() +
geom_col(data = NTJ60.ddPCR %>%
filter(variable != "Mtotal",
Strain == "MUT") %>%
mutate(value = value * 20 / 4 * 40 * 20 / 12 * 60 / 280 * 1000 / 1e7),
aes(x = Locus,
y = value,
fill = variable),
color = "black",
position = "stack") +
facet_wrap(~Rep) +
scale_y_continuous(limits = c(0,10),breaks = 2 * 1:5, labels = comma) +
labs(x = NULL,
y = TeX("\\textbf{vRNA Concentration (copies * $10^7$ / mL)}"),
fill = "Segment") +
scale_fill_manual(values = c("M1.Only" = "steelblue",
"M2.Only" = "lightgoldenrod",
"NS" = "gray40")) +
theme(text=element_text(size=14,face="bold"),
strip.text.x = element_blank(),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),face = "bold",color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
ggsave(g_legend(Legend.7B),
file = "Figures/7BLegend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")
# Figure 7C: Dilution Series ----
NTJ76 = read.csv(file = file.path(Data.Path,"NTJ76_Bayes.csv"), header = TRUE)
NTJ76.Melt = NTJ76 %>%
melt(id = "Dilution") %>%
dplyr::rename(Gate = variable,
Percent = value) %>%
mutate(Marker = substr(Gate,1,2),
Parent = str_c(substr(Gate,4,5),"+"),
Marker_Class = substr(Marker,1,1))
m1 = lm(data = NTJ76.Melt %>% filter(Gate %in% c("HAgM1","HAgM2",
"M1gM2","M2gM1")),
formula = Percent ~ Dilution * Gate)
summary(m1)##
## Call:
## lm(formula = Percent ~ Dilution * Gate, data = NTJ76.Melt %>%
## filter(Gate %in% c("HAgM1", "HAgM2", "M1gM2", "M2gM1")))
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.8200 -2.3133 -0.1917 2.8525 20.5800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.707 2.987 33.717 < 2e-16 ***
## Dilution -10.353 2.439 -4.245 9.03e-05 ***
## GateHAgM2 1.987 4.224 0.470 0.64009
## GateM1gM2 -24.993 4.224 -5.917 2.60e-07 ***
## GateM2gM1 -13.500 4.224 -3.196 0.00237 **
## Dilution:GateHAgM2 1.447 3.449 0.419 0.67661
## Dilution:GateM1gM2 -21.373 3.449 -6.197 9.39e-08 ***
## Dilution:GateM2gM1 -21.340 3.449 -6.188 9.72e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.679 on 52 degrees of freedom
## Multiple R-squared: 0.9504, Adjusted R-squared: 0.9437
## F-statistic: 142.3 on 7 and 52 DF, p-value: < 2.2e-16newdata = matrix(ncol = 4,
nrow = 21 * 4,
data = 0)
colnames(newdata) = c("Dilution","Gate","Parent","Marker_Class")
newdata = data.frame(newdata)
newdata$Dilution = rep(0:20/10,4)
newdata$Gate = rep(c("HAgM1","HAgM2","M1gM2","M2gM1"),each = 21)
newdata$Parent = rep(c("M2","M1"),each = 21)
newdata$Marker_Class = substr(newdata$Marker,1,1)
newdata$Percent = predict(m1,newdata,type="response",se.fit=TRUE)$fit
newdata$se = predict(m1,newdata,type="response",se.fit=TRUE)$se.fit
ggplot() +
geom_jitter(data = NTJ76.Melt %>% filter(Gate %in% c("HAgM1","HAgM2",
"M1gM2","M2gM1")),
aes(x = Dilution,
y = Percent,
color = Gate,
fill = Gate,
shape = Parent),
width = 0.1,
size = 2.5) +
geom_line(data = newdata,
aes(x = Dilution,
y = Percent,
group = Gate,
color = Gate),
lwd = 1.2) +
geom_ribbon(data = newdata,
aes(x = Dilution,
ymin = Percent - 1.96 * se,
ymax = Percent + 1.96 * se,
group = Gate,
fill = Gate),
alpha = 0.3) +
scale_fill_manual(values = c("HAgM1" = "red","HAgM2" = "orange","M1gM2" = "navy","M2gM1" = "skyblue")) +
scale_color_manual(values = c("HAgM1" = "red","HAgM2" = "orange","M1gM2" = "navy","M2gM1" = "skyblue")) +
scale_x_continuous(limits = c(-0.1,2.1)) +
scale_y_continuous(limits = c(-0.5,110),breaks = c(0,25,50,75,100)) +
scale_shape_manual(values = c(24:25)) +
labs(x = TeX("\\textbf{$\\log_1_0$ Stock Dilution}"),
y = TeX("\\textbf{% Protein-Expressing Cells}"),
color = "Protein",
fill = "Protein") +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),face = "bold",color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = c(2,2))
ggsave('Figures/7C_Pan99_M_STOP_Dilution_Series_FACS.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ76.Melt %>% filter(Gate %in% c("HAgM1","HAgM2",
"M1gM2","M2gM1")),
file = file.path(Data.Path,"Fig7C_Data.csv"), row.names = FALSE) # The g stands for "given", as in "HA+ | M1+"
# Figure 7D: WT / M.STOP Ratios ----
Measures = c("NS_ddPCR","M_ddPCR","HA_FACS","PFU","TCID50","GPID50")
NTJ60.Ratios = matrix(ncol = 2,
nrow = length(Measures),
data = 0)
colnames(NTJ60.Ratios) = c("Assay","Ratio")
NTJ60.Ratios = data.frame(NTJ60.Ratios)
NTJ60.Ratios[,"Assay"] = Measures
NTJ60.Ratios[,"Ratio"] = as.numeric(0)
# Ratios 1,2: M and NS RNA copies by ddPCR Titer ----
NTJ60.ddPCR.Sum = NTJ60.ddPCR %>%
na.omit %>%
group_by(variable,Strain) %>%
dplyr::summarise(Copies = mean(value)) %>%
dcast(variable ~ Strain) %>%
filter(variable %in% c("NS","Mtotal")) %>%
mutate(Ratio = WT / MUT)## Using Copies as value column: use value.var to override.NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "NS_ddPCR","Ratio"] = NTJ60.ddPCR.Sum %>% filter(variable == "NS") %>% dplyr::select(Ratio)
NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "M_ddPCR","Ratio"] = NTJ60.ddPCR.Sum %>% filter(variable == "Mtotal") %>% dplyr::select(Ratio)
# Ratio 3: HA FACS Immunotitration ----
NTJ60.FACS = NTJ60 %>%
filter(Assay == "FACS") %>%
mutate(Dilution = -Dilution)
WT.FACS = drm(data = NTJ60 %>% filter(Assay == "FACS",Strain == "WT", Type == "HA") %>% mutate(Dilution = -Dilution),
formula = Readout ~ Dilution,
fct = LL.5()) %>%
ED(50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 1.78977 0.76513M.STOP.FACS = drm(data = NTJ60 %>% filter(Assay == "FACS",Strain == "MUT", Type == "HA") %>% mutate(Dilution = -Dilution),
formula = Readout ~ Dilution,
fct = LL.5()) %>%
ED(50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 1.3770 2.4273NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "HA_FACS","Ratio"] = (10^(WT.FACS - M.STOP.FACS))[,"Estimate"]
# Ratio 4: PFU / mL Titer ----
NTJ60.PFU = NTJ60 %>%
filter(Assay == "PFU")
NTJ60.PFU.Sum = NTJ60.PFU %>%
group_by(Strain) %>%
dplyr::summarise(PFU = mean(Readout))
NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "PFU","Ratio"] =
NTJ60.PFU.Sum %>% filter(Strain == "WT") %>% dplyr::select(PFU) /
NTJ60.PFU.Sum %>% filter(Strain == "MUT") %>% dplyr::select(PFU)
# Ratio 5: TCID50 units / mL Titer ----
NTJ60.TCID50 = NTJ60 %>%
filter(Assay == "TCID50") %>%
mutate(Dilution = -Dilution)
WT.TCID50 = drm(data = NTJ60 %>% filter(Assay == "TCID50",Strain == "WT") %>% mutate(Dilution = -Dilution),
formula = Readout ~ Dilution,
fct = LL.5()) %>%
ED(50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 5.91965 0.35327M.STOP.TCID50 = drm(data = NTJ60 %>% filter(Assay == "TCID50",Strain == "MUT") %>% mutate(Dilution = -Dilution),
formula = Readout ~ Dilution,
fct = LL.5()) %>%
ED(50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 3.838880 0.011153NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "TCID50","Ratio"] = (10^(WT.TCID50 - M.STOP.TCID50))[,"Estimate"]
# WT / M.STOP Ratios for Figure 7D are normalized to NS RNA copy numbers in stocks
# NS was chosen to avoid potentially confounding effects from the M segment mutations
# Note that GPID50 titers are not normalized to relative amounts of NS RNA copies in the analysis because the doses used to generate ID50 curves were in terms of NS RNA copy numbers.
NS_Ratio = (NTJ60.Ratios %>% filter(Assay == "NS_ddPCR") %>% dplyr::select(Ratio) %>% as.numeric)
NTJ60.Ratios = NTJ60.Ratios %>%
mutate(Ratio = as.numeric(Ratio) / NS_Ratio)
# Ratio 6: Guinea Pig ID50 ----
NTJ.GP.ID50 = read.csv(file = file.path(Data.Path,"NTJ_GP_Infected.csv"), header = TRUE) %>%
mutate(Dose_Group = Dose_Group - 1) %>%
group_by(Strain, Dose_Group) %>%
dplyr::summarise(Infected = mean(Infected) * 100)
WT.Model = drm(data = NTJ.GP.ID50 %>% filter(Strain == "WT"),
formula = Infected ~ Dose_Group,
fct = LL.5())
MUT.Model = drm(data = NTJ.GP.ID50 %>% filter(Strain == "MUT"),
formula = Infected ~ Dose_Group,
fct = LL.4())
NTJ60.Ratios[NTJ60.Ratios[,"Assay"] == "GPID50","Ratio"] = as.numeric(10 ^ (WT.Model %>% ED(50) - MUT.Model %>% ED(50)))[1]##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 3.812266 0.076814
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.9010726 0.0045787# Plot Ratios ----
NTJ60.Ratios$Assay = factor(NTJ60.Ratios$Assay,levels = Measures)
ggplot() +
geom_bar(data = NTJ60.Ratios,
stat = "identity",
aes(x = Assay,
y = Ratio %>% as.numeric),
color = "black",
lwd = 1,
fill = "black") +
scale_y_continuous(trans = "log10", limits = c(0.1,1000), breaks = c(1,10,100,1000), labels = comma) +
labs(x = "Method",
y = "Ratio WT / M.STOP") +
coord_cartesian(ylim = c(1,1000)) +
scale_x_discrete(labels = c("NS ddPCR",
"M ddPCR",
"HA FACS",
"PFU",
"TCID ",
"GPID ")) +
theme(text=element_text(size=14,face="bold"),
strip.text.x=element_text(size=rel(1.5),margin=margin(0,0,3,0)),
strip.text.y=element_text(size=rel(1.5),margin=margin(0,0,0,0),angle=0),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=45,hjust = 1,vjust=1,size=rel(1.5),face = "bold",color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "NA")
ggsave('Figures/7D_Pan99_WT_vs_M_STOP_Ratios.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ60.Ratios,
file = file.path(Data.Path,"Fig7D_Data.csv"), row.names = FALSE)
NTJ60.Ratios## Assay Ratio
## 1 NS_ddPCR 1.000000
## 2 M_ddPCR 1.117016
## 3 HA_FACS 1.100739
## 4 PFU 24.021033
## 5 TCID50 51.252994
## 6 GPID50 815.067407# For ddPCR, the limit of detection is 1 cDNA copy / 20 uL reaction, or 0.05 copies / uL
# cDNA copies / uL ddPCR reaction * 20 / 4 * 40 * 20 / 12.8 * 60 / 140 * 1000 = RNA copies / mL sample
RNA_LOD = log10(0.05 *
20 / # 20 uL Reacion
4 * # 4 uL / Reaction
40 * # 1:40 Dilution of cDNA
20 / # 20 uL cDNA / RT reaction
12.8 * # 12.uL RNA / RT reaction
60 / # Total RNA extracted
140 * # uL volume virus used
1000) # uL -> mL
Fig.8 = ggplot() +
labs(y = TeX("\\textbf{Viral Shedding ($\\log_1_0$ RNA copies / mL)}"),
x = "Time Post-Inoculation (Days)") +
geom_hline(yintercept = RNA_LOD,lty = 2) +
scale_y_continuous(limits = c(3.75,8.75),breaks = 4:8) +
theme(text=element_text(size=14,face="bold"),
strip.text.x = element_blank(),
strip.text.y = element_blank(),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
# Figure 8A: WT vs. M.STOP Shedding ----
GP_NS = read.csv(file = file.path(Data.Path,"NTJ69_NTJ73_RNA_Data.csv"), header = TRUE) %>%
mutate(RNA_Shed = log10(RNA_Conc *
20 / # 20 uL Reacion
4 * # 4 uL / Reaction
40 * # 1:40 Dilution
20 /
12.8 * # RNA / cDNA reaction
60 / # Total RNA extracted
140 * # volume virus used
1000), # uL -> mL
Virus = recode(Strain,"MUT" = "M.STOP") %>% factor(levels = c("WT","M.STOP")))
Fig.8 +
geom_line(data = GP_NS,
aes(x = Day,
y = RNA_Shed,
color = Virus,
group = GP),
lwd = 1.1) +
labs(y = TeX("\\textbf{Viral Shedding ($\\log_1_0$ RNA copies / mL)}"),
x = "Time Post-Inoculation (Days)") +
scale_x_continuous(breaks = c(1,2,3,5,7)) +
scale_color_manual(values = c("WT" = "black",
"M.STOP" = "mediumseagreen")) +
geom_hline(yintercept = RNA_LOD,lty = 2) +
theme(legend.position = "NA")
ggsave('Figures/8A_Pan99_WT_M_STOP_Shedding_Comparison.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(GP_NS,
file = file.path(Data.Path,"Fig8A_Data.csv"), row.names = FALSE)
Legend.8A = Fig.8 +
geom_line(data = GP_NS,
aes(x = Day,
y = RNA_Shed,
color = Virus,
group = GP),
lwd = 1.1) +
labs(y = TeX("\\textbf{Viral Shedding ($\\log_1_0$ RNA copies / mL)}"),
x = "Time Post-Inoculation (Days)") +
scale_x_continuous(breaks = c(1,2,3,5,7)) +
scale_color_manual(values = c("WT" = "black",
"M.STOP" = "mediumseagreen")) +
geom_hline(yintercept = RNA_LOD,lty = 2)
Fig.8A.Legend = arrangeGrob(g_legend(Legend.8A))
ggsave(Fig.8A.Legend,
file = "Figures/8A_Legend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")
# Figure 8B: WT Dose Comparison----
NTJ77 = read.csv(file = file.path(Data.Path,"NTJ77_Data.csv"), header = TRUE)
Neg_RNA = NTJ77 %>% filter(Dose == "Neg") %>% dplyr::select(RNA_Conc) %>% as.numeric
NTJ77 = NTJ77 %>%
mutate(RNA_Shed = log10(pmax(Neg_RNA, RNA_Conc) *
20 / # 20 uL Reacion
4 * # 4 uL / Reaction
40 * # 1:40 Dilution of cDNA
20 / # 20 uL cDNA / RT reaction
12.8 * # RNA / RT reaction
60 / # Total RNA extracted
140 * # volume virus used
1000)) # uL -> mL
Fig.8 +
geom_line(data = NTJ77 %>% filter(Dose != "Neg"),
aes(x = Day,
y = RNA_Shed,
group = GP,
lty = Dose),
lwd = 1.2) +
labs(linetype = "WT Dose") +
scale_x_continuous(breaks = c(1,2,3,5,7)) +
theme(legend.position = "NA")
ggsave('Figures/8B_Pan99_WT_Shedding_High_vs_Low_Dose.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(NTJ77 %>% filter(Dose != "Neg"),
file = file.path(Data.Path,"Fig8B_Data.csv"), row.names = FALSE)
Legend.8B = Fig.8 +
geom_line(data = NTJ77 %>% filter(Dose != "Neg"),
aes(x = Day,
y = RNA_Shed,
group = GP,
lty = Dose),
lwd = 1.2) +
labs(linetype = "WT Dose") +
scale_x_continuous(breaks = c(1,2,3,5,7))
Fig.8B.Legend = arrangeGrob(g_legend(Legend.8B))
ggsave(Fig.8B.Legend,
file = "Figures/8B_Legend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")
# Figure 8C: Transmission ----
NTJ74.Short = read.csv(file = file.path(Data.Path,"NTJ74_ddPCR_Data.csv"), header = TRUE)
NTJ74.Long = NTJ74.Short %>% melt(id = c("GP","Pair","Cage","Status","Strain")) %>%
dplyr::rename(Day = variable,
RNA = value) %>%
mutate(Day = as.factor(Day) %>% recode_factor("X2" = 2, "X4" = 4, "X6" = 6, "X8" = 8) %>% paste %>% as.numeric,
Strain = str_replace(Strain,"MUT","M.STOP"),
Strain = factor(Strain,levels = c("WT","M.STOP")),
RNA = log10(RNA *
20 / # 20 uL Reacion
4 * # 4 uL cDNA/ Reaction
40 * # 1:40 Dilution
20 /
12.8 * # RNA / cDNA reaction
60 / # Total RNA extracted
140 * # volume virus used
1000), # uL -> mL)
Virus = Strain)
Fig.8 +
geom_line(data = NTJ74.Long %>% filter(Strain == "WT"),
aes(x = Day,
y = RNA,
linetype = Status,
color = Virus,
group = GP),
lwd = 1.1) +
facet_grid(~ Cage) +
scale_x_continuous(breaks = c(2,4,6,8)) +
scale_color_manual(values = c("WT" = "black",
"M.STOP" = "mediumseagreen")) +
geom_hline(yintercept = RNA_LOD,lty = 2) +
theme(legend.position = "NA")
ggsave('Figures/8C_Pan99_WT_vs_M_STOP_Transmission_RNA.pdf',
width = 10,
height = 5,
unit = "in")
write.csv(NTJ74.Long %>% filter(Strain == "WT"),
file = file.path(Data.Path,"Fig8C_Data.csv"), row.names = FALSE)
Fig.8 +
geom_line(data = NTJ74.Long %>% filter(Strain == "M.STOP"),
aes(x = Day,
y = RNA,
linetype = Status,
color = Virus,
group = GP),
lwd = 1.1) +
facet_grid(~ Cage) +
scale_x_continuous(breaks = c(2,4,6,8)) +
scale_color_manual(values = c("WT" = "black",
"M.STOP" = "mediumseagreen")) +
geom_hline(yintercept = RNA_LOD,lty = 2) +
theme(legend.position = "NA")
ggsave('Figures/8D_Pan99_WT_vs_M_STOP_Transmission_RNA.pdf',
width = 10,
height = 5,
unit = "in")
write.csv(NTJ74.Long %>% filter(Strain == "M.STOP"),
file = file.path(Data.Path,"Fig8D_Data.csv"), row.names = FALSE)
Legend.Fig.8C = Fig.8 +
geom_line(data = NTJ74.Long,
aes(x = Day,
y = RNA,
linetype = Status,
color = Virus,
group = GP),
lwd = 1.1) +
facet_grid(~ Cage) +
scale_x_continuous(breaks = c(2,4,6,8)) +
scale_color_manual(values = c("WT" = "black",
"M.STOP" = "mediumseagreen")) +
geom_hline(yintercept = RNA_LOD,lty = 2)
Fig.8C.Legend = arrangeGrob(g_legend(Legend.Fig.8C))
ggsave(Fig.8C.Legend,
file = "Figures/8Legend.pdf",
dpi = 300,
width = 4,
height = 3,
unit = "in")
NTJ69.ddPCR = read.csv(file = file.path(Proj.Home,"Data","NTJ69_ddPCR_Data.csv")) %>%
mutate(M = M1 + M2,
M1_Frac = M1 / M,
M2_Frac = M2 / M,
M_NS = M / NS,
M = log10(M * 20 / 4 * 40 * 20 / 12.8 * 60 / 140 * 1000),
M1 = log10(M1 * 20 / 4 * 40 * 20 / 12.8 * 60 / 140 * 1000),
M2 = log10(M2 * 20 / 4 * 40 * 20 / 12.8 * 60 / 140 * 1000),
NS = log10(NS * 20 / 4 * 40 * 20 / 12.8 * 60 / 140 * 1000))
NTJ60.ddPCR = read.csv(file = file.path(Proj.Home,"Data","NTJ60_Data.csv"), header = TRUE) %>% filter(Assay == "ddPCR",
Type %in% c("Mwt","M2stop","NS")) %>%
dcast(Strain + Rep ~ Type) %>%
mutate(Mtotal = M2stop + Mwt) %>%
melt(id = c("Strain","Rep")) %>%
na.omit %>%
mutate(Locus = str_replace_all(variable,c("Mwt" = "M", "M2stop" = "M")),
variable = str_replace_all(variable,c("Mwt" = "M2.Only", "M2stop" = "M1.Only")))## Using Readout as value column: use value.var to override.M1 = NTJ60.ddPCR %>% filter(variable == "M1.Only", Strain == "MUT") %>% select(value)
M2 = NTJ60.ddPCR %>% filter(variable == "M2.Only", Strain == "MUT") %>% select(value)
NS = NTJ60.ddPCR %>% filter(variable == "NS", Strain == "MUT") %>% select(value)
Stock_M1_Frac = mean((M1 /(M1 + M2))$value)
Stock_M2_Frac = mean((M2 /(M1 + M2))$value)
Stock_M_NS = mean(((M1 + M2)/(NS))$value)
Stock.ddPCR = NTJ69.ddPCR %>% filter(Day == 1) %>%
mutate(Day = 0,
M1 = 0,
M2 = 0,
NS = 0,
M = 0,
M1_Frac = Stock_M1_Frac,
M2_Frac = Stock_M2_Frac,
M_NS = Stock_M_NS)
Rescue.ddPCR = Stock.ddPCR %>%
mutate(Day = -1,
M1_Frac = 0.5,
M2_Frac = 0.5,
M_NS = 1)
MSTOP.ddPCR = rbind(Rescue.ddPCR,Stock.ddPCR,NTJ69.ddPCR)
ggplot() +
geom_segment(aes(x = -0.5, y = .9, xend = -0, yend = Stock_M1_Frac),lty = 2) +
geom_segment(aes(x = -0.5, y = .9, xend = -0, yend = Stock_M2_Frac),lty = 2) +
geom_segment(aes(x = -1, y = 0.05, xend = -1, yend = 0.5),lty = 2) +
geom_text(aes(x = -1,y = 0,label = "Rescue",fontface = "bold", size = 14)) +
geom_text(aes(x = -0.5,y = 0.95,label = "Inoculum",fontface = "bold", size = 14)) +
geom_line(data = MSTOP.ddPCR,
aes(x = Day,
y = M1_Frac,
group = GP),
lwd = 1.2,
color = "steelblue") +
geom_line(data = MSTOP.ddPCR,
aes(x = Day,
y = M2_Frac,
group = GP),
lwd = 1.2,
color = "lightgoldenrod") +
geom_point(data = Stock.ddPCR %>% melt(id = c("GP","Day")) %>% filter(variable %in% c("M1_Frac","M2_Frac")),
aes(x = Day,
y = value,
fill = variable),
pch = 21,
stroke = 2,
size = 5) +
geom_point(data = Rescue.ddPCR,
aes(x = Day,
y = M1_Frac),
fill = "mediumseagreen",
pch = 21,
stroke = 2,
size = 5) +
labs(x = "Time (Days Post-Inoculation)",
y = "Proportion of Total M Segments") +
scale_fill_manual(values = c("M1_Frac" = "steelblue", "M2_Frac" = "lightgoldenrod")) +
scale_x_continuous(breaks = c(0,1,2,3),
limits = c(-1.25,3)) +
scale_y_continuous(limits = c(0,1)) +
theme(text=element_text(size=14,face="bold"),
strip.text.x = element_blank(),
strip.text.y = element_blank(),
strip.background = element_blank(),
plot.title = element_text(size=rel(1.5)),
axis.text.x=element_text(angle=0,vjust=0,size=rel(1.5),color = "black"),
axis.text.y=element_text(size=rel(1.5),color = "black"),
axis.line.x = element_line(size=0.5),
axis.line.y = element_line(size=0.5),
axis.ticks.x = element_line(size=0.5, color = "black"),
axis.ticks.y = element_line(size = 0.5, color = "black"),
axis.title.y = element_text(size=rel(1.2),color = "black"),
axis.title.x = element_text(size=rel(1.2)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = c(2,2))
ggsave('Figures/Supp5B_M1_M2_Ratio.pdf',
width = 5,
height = 5,
unit = "in")
write.csv(MSTOP.ddPCR,
file = file.path(Data.Path,"Supp5B_Data.csv"), row.names = FALSE)